Dear DuMux Community, I hope this email finds you well and safe.
I am writing to kindly ask for your support in the following:
a) I am working on porousmedium flow model (injection problem)
b)I have a domain of (10m x 10m) and it includes an aquifer, cement layer, and
overburden (as you can see in the attached image1)
c) all layers have a constant porosity and permeability unless the cement layer
d) the cement layer is defined as following
bool isInCement_(const GlobalPosition &globalPos) const
{
return globalPos[dimWorld-2] > 5 + eps_ && globalPos[dimWorld-2] <
5.2+ eps_ && globalPos[dimWorld-1] > 0 + eps_ && globalPos[dimWorld-1] < 7+
eps_; }
e) the cement layer has a length of 20 cm and the porosity changes as function
of globalPos[2] and time
f) the cement layer permeability is dependent on porosity
g) I am using cctpfa for discretization
h) the grid as following:
Positions0 = 0 2 5 5.2 10
Positions1 = 0 10
Cells0 = 100 100 100 100
Cells1 = 200
Grading0 = 1.0 1.0 1.0 1.0
Grading1 = 1.0
I) the resolution in the cement layer (x-direction) should be small (e.g.
0.001).
Issue:
a) at the beginning of the simulation I get correct porosity values of the
cement but after few minutes I start to get very large porosity values (either
positive or negative) and I am not able to understand the reason behind that.
note: when I change the position of the cement layer from the middle of the
domain (5-5.2) to the left of the domain (0-0.2) as you can see below, I got
the right porosity values of the cement layer (image 2).
bool isInCement_(const GlobalPosition &globalPos) const
{
return globalPos[dimWorld-2] > 0 + eps_ && globalPos[dimWorld-2] <
0.2+ eps_ && globalPos[dimWorld-1] > 0 + eps_ && globalPos[dimWorld-1] < 7+
eps_; }
Positions0 = 0 0.2 5 10
Positions1 = 0 10
Cells0 = 100 100 100 100
Cells1 = 200
Grading0 = 1.0 1.0 1.0 1.0
Grading1 = 1.0
I attached the files for things to be clear.
Looking forward to your help.
Best Regards,
Mohammad Hodroj
params.input
Description: params.input
// -*- 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 two-phase porousmediumflow properties file for exercise-basic */ #ifndef DUMUX_EX_BASIC_PROPERTIES_2P_HH #define DUMUX_EX_BASIC_PROPERTIES_2P_HH #include <dune/grid/yaspgrid.hh> #include <dumux/discretization/cctpfa.hh> #include <dumux/discretization/elementsolution.hh> #include <dumux/discretization/method.hh> #include <dumux/discretization/box.hh> #include <dumux/discretization/elementsolution.hh> #include <dumux/porousmediumflow/2p/model.hh> #include <dumux/material/fluidsystems/h2och4.hh> #include <dumux/discretization/box/gridvolumevariables.hh> #include <dumux/discretization/box/fvgridgeometry.hh> #include "injection2pproblem.hh" #include "spatialparams.hh" namespace Dumux::Properties { // define the TypeTag for this problem with a cell-centered two-point flux approximation spatial discretization. // Create new type tags namespace TTag { struct Injection2p { using InheritsFrom = std::tuple<TwoP>; }; struct Injection2pCC { using InheritsFrom = std::tuple<Injection2p, CCTpfaModel>; }; //struct Injection2pBox { using InheritsFrom = std::tuple<Injection2p, BoxModel>; }; } // end namespace TTag // Set the grid type template<class TypeTag> //struct Grid<TypeTag, TTag::Injection2p> { using type = Dune::YaspGrid<2>; }; struct Grid<TypeTag, TTag::Injection2p> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<double, 2>>; }; // Set the problem property template<class TypeTag> struct Problem<TypeTag, TTag::Injection2p> { using type = Injection2PProblem<TypeTag>; }; // Set the spatial parameters template<class TypeTag> struct SpatialParams<TypeTag, TTag::Injection2p> { private: using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; using Scalar = GetPropType<TypeTag, Properties::Scalar>; public: using type = InjectionSpatialParams<FVGridGeometry, Scalar>; }; // Set fluid configuration template<class TypeTag> struct FluidSystem<TypeTag, TTag::Injection2p> { using type = FluidSystems::H2OCH4< GetPropType<TypeTag, Properties::Scalar>, FluidSystems::H2OCH4DefaultPolicy</*fastButSimplifiedRelations=*/ true> >; }; } // end namespace Dumux::Properties #endif
// -*- 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 Definition of the spatial parameters for the injection problem * which uses the isothermal two-phase two-component * fully implicit model. */ #ifndef DUMUX_EX_BASIC_SPATIAL_PARAMS_HH #define DUMUX_EX_BASIC_SPATIAL_PARAMS_HH #include <dumux/material/spatialparams/fv.hh> #include <dumux/material/fluidmatrixinteractions/2p/brookscorey.hh> #include <cmath> #include <dumux/io/gnuplotinterface.hh> #include <dumux/io/plotpckrsw.hh> namespace Dumux { /*! * \ingroup TwoPTwoCModel * \brief Definition of the spatial parameters for the injection problem * which uses the isothermal two-phase two-component * fully implicit model. */ template<class FVGridGeometry, class Scalar> class InjectionSpatialParams : public FVSpatialParams<FVGridGeometry, Scalar, InjectionSpatialParams<FVGridGeometry, Scalar>> { using ThisType = InjectionSpatialParams<FVGridGeometry, Scalar>; using ParentType = FVSpatialParams<FVGridGeometry, Scalar, ThisType>; using GridView = typename FVGridGeometry::GridView; using FVElementGeometry = typename FVGridGeometry::LocalView; // get the dimensions of the simulation domain from GridView static const int dimWorld = GridView::dimensionworld; using Element = typename GridView::template Codim<0>::Entity; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; using PcKrSwCurve = FluidMatrix::BrooksCoreyDefault<Scalar>; public: // export permeability type using PermeabilityType = Scalar; /*! * \brief The constructor * * \param fvGridGeometry The finite volume grid geometry */ InjectionSpatialParams(std::shared_ptr<const FVGridGeometry>& fvGridGeometry) : ParentType(fvGridGeometry) , cementPcKrSwCurve_("SpatialParams.Cement") , overburdenPcKrSwCurve_("SpatialParams.Overburden") , aquiferPcKrSwCurve_("SpatialParams.Aquifer") { // intrinsic permeabilities overburdenK_ = getParam<Scalar>("SpatialParams.PermeabilityOverburden"); aquiferK_ = getParam<Scalar>("SpatialParams.PermeabilityAquifer"); // porosities overburdenPorosity_ = 0.01; aquiferPorosity_ = 0.33; } /* * \brief Define the intrinsic permeability \f$\mathrm{[m^2]}\f$. * * \param globalPos The global position */ template<class GlobalPosition> PermeabilityType permeabilityAtPos (const GlobalPosition& globalPos)const // here, either aquitard or aquifer permeability are returned, depending on the global position { Scalar cementPor_; Scalar cementK_; const Scalar dc_Porosity = 0.3; Scalar InitialCementK_ = 8.645e-18; if (isInCement_(globalPos)) { cementPor_ = porosityAtPos( globalPos); cementK_ = InitialCementK_ *(std::pow(((1-dc_Porosity)/(1- cementPor_)),2) *(std::pow(( cementPor_/dc_Porosity), 3))); std::cout<< " Porosity:"<< cementPor_ <<" Permeability:"<< cementK_ << std::endl; return cementK_; } else if (isInAquifer_(globalPos)) return aquiferK_; else return overburdenK_; } /* * \brief Define the porosity \f$\mathrm{[-]}\f$. * * \param globalPos The global position */ Scalar porosityAtPos(const GlobalPosition& globalPos) const { // here, either aquitard or aquifer porosity are returned, depending on the global position if (isInCement_(globalPos)) { // constants to define const Scalar x_min_cement = 5; const Scalar x_max_cement = 5.2; const Scalar dc_Porosity = 0.3; // Scalar fc_porosity; Scalar slope_porosity; Scalar cementPorosity_; Scalar dc1 = 0.0054*(std::sqrt(time_/31536000.0))+0.001; Scalar fc1 = dc1/1.375; Scalar fc = x_min_cement + fc1; Scalar dc = x_min_cement + dc1; using std::pow; if (((globalPos[dimWorld-2]) >= x_min_cement) && ((globalPos[dimWorld-2]) <= (std::min(fc, x_max_cement)))){ Scalar a = -0.0005 * std::pow((time_/31536000.0), 2)+ 1.3776 * (time_/31536000.0) - 987.37; Scalar b = -7e-5 * std::pow((time_/31536000.0), 2)+ 0.0228 * (time_/31536000.0) +87.137; Scalar c = 8e-6 * std::pow((time_/31536000.0), 2) - 0.0102 * (time_/31536000.0) -1.9432; Scalar d = -9e-8 * std::pow((time_/31536000.0), 2) + 0.0003 * (time_/31536000.0) + 0.3804; cementPorosity_ = a*std::pow(globalPos[dimWorld-2],3) + b*std::pow(globalPos[dimWorld-2],2) + c*globalPos[dimWorld-2] + d; } else if ((globalPos[dimWorld-2]) >= fc && ((globalPos[dimWorld-2]) <= std::min(dc, x_max_cement))){ Scalar a = -0.0005 * std::pow((time_/31536000.0), 2)+ 1.3776 * (time_/31536000.0) - 987.37; Scalar b = -7e-5 * std::pow((time_/31536000.0), 2)+ 0.0228 * (time_/31536000.0) +87.137; Scalar c = 8e-6 * std::pow((time_/31536000.0), 2) - 0.0102 * (time_/31536000.0) -1.9432; Scalar d = -9e-8 * std::pow((time_/31536000.0), 2) + 0.0003 * (time_/31536000.0) + 0.3804; this->fc_porosity = a*std::pow(fc,3) + b*std::pow(fc,2) + c*fc + d; Scalar slope_porosity = (dc_Porosity - this->fc_porosity)/(dc - fc); cementPorosity_ = std::abs(slope_porosity *(globalPos[dimWorld-2]-fc)+fc_porosity); } else { cementPorosity_ = dc_Porosity; } return cementPorosity_; } else if (isInAquifer_(globalPos)) return aquiferPorosity_; else return overburdenPorosity_; } /*! * \brief Returns the fluid-matrix interaction law at a given location * \param globalPos The global coordinates for the given location */ auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const { if (isInCement_(globalPos)) return makeFluidMatrixInteraction(cementPcKrSwCurve_); else if (isInAquifer_(globalPos)) return makeFluidMatrixInteraction(aquiferPcKrSwCurve_); else return makeFluidMatrixInteraction(overburdenPcKrSwCurve_); } /*! * \brief Function for defining which phase is to be considered as the wetting phase. * * \return the wetting phase index * \param globalPos The position of the center of the element */ template<class FluidSystem> int wettingPhaseAtPos(const GlobalPosition& globalPos) const { return FluidSystem::H2OIdx; } //! set the time for the time dependent boundary conditions (called from main) void setTime(Scalar time) { time_ = time; /*std::cout<<"Helloz from time spatial\n";*/ } private: Scalar time_; mutable Scalar fc_porosity; //mutable Scalar cementK_; static constexpr Scalar eps_ = 1e-6; // provides a convenient way distinguishing whether a given location is inside the aquitard bool isInCement_(const GlobalPosition &globalPos) const { // globalPos[dimWorld-1] is the y direction for 2D grids or the z direction for 3D grids return globalPos[dimWorld-2] > 5 + eps_ && globalPos[dimWorld-2] < 5.2+ eps_ && globalPos[dimWorld-1] > 0 + eps_ && globalPos[dimWorld-1] < 7+ eps_; } bool isInAquifer_(const GlobalPosition &globalPos) const { // globalPos[dimWorld-1] is the y direction for 2D grids or the z direction for 3D grids return globalPos[dimWorld-2] > 0+ eps_ && globalPos[dimWorld-2] < 10- eps_ && globalPos[dimWorld-1] > 7 + eps_ && globalPos[dimWorld-1] < 10 + eps_; } Scalar overburdenK_; //Scalar cementzK_; Scalar aquiferK_; Scalar cementPorosity_; Scalar overburdenPorosity_; Scalar aquiferPorosity_; const PcKrSwCurve cementPcKrSwCurve_; const PcKrSwCurve overburdenPcKrSwCurve_; const PcKrSwCurve aquiferPcKrSwCurve_; }; } // end namespace Dumux #endif
// -*- 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 main file for the two-phase porousmediumflow problem of exercise-basic */ #include <config.h> #include <ctime> #include <iostream> #include <dune/common/parallel/mpihelper.hh> #include <dune/common/timer.hh> #include <dune/grid/io/file/dgfparser/dgfexception.hh> #include <dune/grid/io/file/vtk.hh> #include <dumux/common/properties.hh> #include <dumux/common/parameters.hh> #include <dumux/common/dumuxmessage.hh> #include <dumux/common/defaultusagemessage.hh> #include <dumux/linear/amgbackend.hh> #include <dumux/linear/linearsolvertraits.hh> #include <dumux/nonlinear/newtonsolver.hh> #include <dumux/assembly/fvassembler.hh> #include <dumux/assembly/diffmethod.hh> #include <dumux/discretization/method.hh> #include <dumux/io/vtkoutputmodule.hh> #include <dumux/io/grid/gridmanager.hh> // The properties file, where the compile time options are defined #include "properties2p.hh" //////////////////////// // the main function //////////////////////// int main(int argc, char** argv) { using namespace Dumux; // define the type tag for this problem using TypeTag = Properties::TTag::Injection2pCC; // initialize MPI, finalize is done automatically on exit const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); // print dumux start message if (mpiHelper.rank() == 0) DumuxMessage::print(/*firstCall=*/true); // parse command line arguments and input file Parameters::init(argc, argv); // try to create a grid (from the given grid file or the input file) GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager; gridManager.init(); //////////////////////////////////////////////////////////// // run instationary non-linear problem on this grid //////////////////////////////////////////////////////////// // we compute on the leaf grid view const auto& leafGridView = gridManager.grid().leafGridView(); // create the finite volume grid geometry using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView); fvGridGeometry->update(); // the problem (initial and boundary conditions) using Problem = GetPropType<TypeTag, Properties::Problem>; auto problem = std::make_shared<Problem>(fvGridGeometry); // the solution vector using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; //SolutionVector x; SolutionVector x(fvGridGeometry->numDofs()); problem->applyInitialSolution(x); auto xOld = x; //problem->spatialParams().permeabilityAtPos (*fvGridGeometry, x); // the grid variables using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry); gridVariables->init(x); // get some time loop parameters using Scalar = GetPropType<TypeTag, Properties::Scalar>; const auto tEnd = getParam<Scalar>("TimeLoop.TEnd"); const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize"); auto dt = getParam<Scalar>("TimeLoop.DtInitial"); // intialize the vtk output module using IOFields = GetPropType<TypeTag, Properties::IOFields>; // use non-conforming output for the test with interface solver const auto ncOutput = getParam<bool>("Problem.UseNonConformingOutput", false); VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name(), "", ncOutput ? Dune::VTK::nonconforming : Dune::VTK::conforming); using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>; vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables)); GetPropType<TypeTag, Properties::IOFields>::initOutputModule(vtkWriter); //!< Add model specific output fields vtkWriter.addField(problem->getpermeability(), "Permeability"); problem->updateVtkOutput(x); vtkWriter.write(0.0); // instantiate time loop auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd); timeLoop->setMaxTimeStepSize(maxDt); // the assembler with time loop for instationary problem using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>; auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop, xOld); // the linear solver using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<FVGridGeometry>>; auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper()); // the non-linear solver using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>; NewtonSolver nonLinearSolver(assembler, linearSolver); // time loop timeLoop->start(); while (!timeLoop->finished()) { //std::cout<<"\n\nMAIN\n\n"; // set previous solution for storage evaluations assembler->setPreviousSolution(xOld); //set time in problem (is used in time-dependent Neumann boundary condition) problem->setTime(timeLoop->time()+timeLoop->timeStepSize()); // solve the non-linear system with time step control nonLinearSolver.solve(x, *timeLoop); // make the new solution the old solution xOld = x; gridVariables->advanceTimeStep(); // advance to the time loop to the next step timeLoop->advanceTimeStep(); // update the output fields before write problem->updateVtkOutput(x); // report statistics of this time step timeLoop->reportTimeStep(); // set new dt as suggested by the newton solver timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize())); problem->updateVtkOutput(x); // output to vtk vtkWriter.write(timeLoop->time()); } timeLoop->finalize(leafGridView.comm()); //////////////////////////////////////////////////////////// // finalize, print dumux message to say goodbye //////////////////////////////////////////////////////////// // print dumux end message if (mpiHelper.rank() == 0) { Parameters::print(); DumuxMessage::print(/*firstCall=*/false); } return 0; } // end main
_______________________________________________ DuMux mailing list [email protected] https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
