Hello, by changing the initially 2p problem for geothermal simulation (made by Samuel Scherrer) into a 3p3c problem and running it, the following error messages occured: "Solve: M deltax^k = rCould not build any aggregates. Probably no connected nodes." For more detailed output please look into output.txt. We tried to change the solvers and to simplify the problem, but nothing worked, we do not know, what the problem could be.
Does anyone know about something about this message: "Could not build any aggregates. Probably no connected nodes." ? Best wishes, Kai Wendel
Let's get the cow off the ice.
Reading parameters from file params.input.
Reading 3d Gmsh grid...
version 2.2 Gmsh file detected
file contains 5097 nodes
file contains 26744 elements
number of real vertices = 5097
number of boundary elements = 7154
number of elements = 19590
Build neighbors took 0.0291691 sec.
NotStrongCompatibleMacroFaces InnerFaces TotalFaces Maximum/Vertex
Minimum/Vertex
15961 35603 42757 36 0
Marking longest edge for initial refinement...
WARNING (ignored): Could not open file 'alugrid.cfg', using default values 0 <
[balance] < 1.2, partitioning method 'ALUGRID_SpaceFillingCurve(9)'.
You are using DUNE-ALUGrid, please don't forget to cite the paper:
Alkaemper, Dedner, Kloefkorn, Nolte. The DUNE-ALUGrid Module, 2016.
Created parallel ALUGrid<3,3,simplex,conforming> from input stream.
Writing output for problem "heat_source". Took 0.534159 seconds.
Enabled periodic check points every 86400 seconds with the next check point at
86400 seconds.
Newton solver configured with the following options and parameters:
-- Newton.EnableShiftCriterion = true (relative shift convergence criterion)
-- Newton.MaxRelativeShift = 0.0001
-- Newton.MinSteps = 2
-- Newton.MaxSteps = 6
-- Newton.TargetSteps = 10
-- Newton.RetryTimeStepReductionFactor = 0.5
-- Newton.MaxTimeStepDivisions = 10
Solve: M deltax^k = rCould not build any aggregates. Probably no connected
nodes.
UMFPACK V5.7.6 (May 4, 2016): WARNING: matrix is singular
Newton: Caught exception: "SolverAbort
[step:/home/kaiw/DUMUX/dune-istl/dune/istl/solver.hh:429]:
Dune::BiCGSTABSolver: defect=-nan is infinite or NaN"
Newton solver did not converge with dt = 10 seconds. Retrying with time step of
5 seconds
Solve: M deltax^k = rCould not build any aggregates. Probably no connected
nodes.
// -*- 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 3 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 * \ingroup VegasGeothermal * \brief 1D problem setup for heat storage simulation in the unsaturated zone * * This is only to generate a stationary initial condition and boundary condition. */ #ifndef DUMUX_VEGAS_GEOTHERMAL_2D_HEATSOURCE_PROBLEM_HH #define DUMUX_VEGAS_GEOTHERMAL_2D_HEATSOURCE_PROBLEM_HH #include <typeinfo> #include <dune/alugrid/grid.hh> #include <dumux/common/variablelengthspline_.hh> #include <dumux/common/intersectionmapper.hh> #include <dumux/discretization/cctpfa.hh> #include <dumux/discretization/box.hh> #include <dumux/discretization/method.hh> #include <dumux/discretization/evalsolution.hh> #include <dumux/discretization/evalgradients.hh> #include <dumux/discretization/box/scvftoscvboundarytypes.hh> // Added CCMPFA.hh! #include <dumux/discretization/ccmpfa.hh> #include <dumux/material/binarycoefficients/h2o_air.hh> #include <dumux/material/components/tabulatedcomponent.hh> #include <dumux/material/fluidsystems/h2oairmesitylene.hh> #include <dumux/material/solidstates/compositionalsolidstate.hh> #include <dumux/material/solidsystems/compositionalsolidphase.hh> #include <dumux/material/components/constant.hh> #include <dumux/porousmediumflow/3p3c/model.hh> #include <dumux/porousmediumflow/problem.hh> #include "spatialparams.hh" #include "utils.hh" #include "io.hh" #include "mysimpleh2o.hh" #include "heatsource.hh" namespace Dumux { template <class TypeTag> class HeatSourceProblem; namespace Properties { // Create new type tags namespace TTag { struct HeatSource { using InheritsFrom = std::tuple<ThreePThreeCNI>; }; struct HeatSourceBox { using InheritsFrom = std::tuple<HeatSource, BoxModel>; }; //KAIW Added struct HeatSourceCCMpfaTypeTag {using InheritsFrom = std::tuple<CCMpfaModel ,HeatSource>; }; } // end namespace TTag template<class TypeTag> struct EnableFVGridGeometryCache<TypeTag, TTag::HeatSource> { static constexpr bool value = true; }; template<class TypeTag> struct EnableGridFluxVariablesCache<TypeTag, TTag::HeatSource> { static constexpr bool value = true; }; template<class TypeTag> struct EnableGridVolumeVariablesCache<TypeTag, TTag::HeatSource> { static constexpr bool value = true; }; // Set the grid type template<class TypeTag> struct Grid<TypeTag, TTag::HeatSource>{ using type = Dune::ALUGrid<3, 3, Dune::simplex, Dune::conforming>; }; /* Since boundary flags/boundary domain markers work differently than normal for * ALUGrid, we need to use a BoundarySegmentIndexFlag instead of BoundaryFlag * for the SCVFs. * Therefore, we first have to define custom grid geometry traits and then set a * non-standard FVGridGeometry property. * This depends on the chosen discretization method, in this case Box. */ template<class GridView> struct MyGridGeometryTraits : public BoxDefaultGridGeometryTraits<GridView> { struct MyScvfTraits : public BoxDefaultScvfGeometryTraits<GridView> { using BoundaryFlag = BoundarySegmentIndexFlag; }; using SubControlVolumeFace = BoxSubControlVolumeFace<GridView, MyScvfTraits>; }; template<class TypeTag> struct FVGridGeometry<TypeTag, TTag::HeatSourceBox> { private: static constexpr bool enableCache = getPropValue<TypeTag, Properties::EnableFVGridGeometryCache>(); using GridView = GetPropType<TypeTag, Properties::GridView>; using Scalar = GetPropType<TypeTag, Properties::Scalar>; public: using type = BoxFVGridGeometry<Scalar, GridView, enableCache, MyGridGeometryTraits<GridView>>; }; // Set the problem property template<class TypeTag> struct Problem<TypeTag, TTag::HeatSource> { using type = HeatSourceProblem<TypeTag>; }; // TODO: Change the names or not, look how it is better // Set the wetting phase template<class TypeTag> struct FluidSystem<TypeTag, TTag::HeatSource> { using Scalar = GetPropType<TypeTag, Properties::Scalar>; using type = FluidSystems::H2OAirMesitylene<Scalar,/*H2Otype=*/DumuxVegasGeothermal::MySimpleH2O<Scalar> >; // seems to work :-) }; // Taken from kuvette! template<class TypeTag> struct SolidSystem<TypeTag, TTag::HeatSource> { using Scalar = GetPropType<TypeTag, Properties::Scalar>; using ComponentOne = Dumux::Components::Constant<1, Scalar>; using ComponentTwo = Dumux::Components::Constant<2, Scalar>; static constexpr int numInertComponents = 2; using type = SolidSystems::CompositionalSolidPhase<Scalar, ComponentOne, ComponentTwo, numInertComponents>; }; template<class TypeTag> struct SolidState<TypeTag, TTag::HeatSource> { private: using Scalar = GetPropType<TypeTag, Properties::Scalar>; using SolidSystem = GetPropType<TypeTag, Properties::SolidSystem>; public: using type = CompositionalSolidState<Scalar, SolidSystem>; }; // Define whether mole(true) or mass (false) fractions are used template<class TypeTag> struct UseMoles<TypeTag, TTag::HeatSource> { static constexpr bool value = true; }; // Set the spatial parameters template<class TypeTag> struct SpatialParams<TypeTag, TTag::HeatSource> { using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>; using Scalar = GetPropType<TypeTag, Properties::Scalar>; using type = DumuxVegasGeothermal::HorizontallyLayeredSpatialParams<FVGridGeometry, Scalar>; }; } // end namespace Properties /*! * \ingroup 2dGeothermal * \brief Setup for a 2d model for unsaturated geothermal storage * * */ template <class TypeTag> class HeatSourceProblem : public PorousMediumFlowProblem<TypeTag> { using ParentType = PorousMediumFlowProblem<TypeTag>; using Scalar = GetPropType<TypeTag, Properties::Scalar>; using GridView = GetPropType<TypeTag, Properties::GridView>; using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; using Indices = typename ModelTraits::Indices; static constexpr int dimWorld = GridView::dimensionworld; // primary variable indices static constexpr int pressureIdx = Indices::pressureIdx; // gas phase pressure //--------// new copied from NAPL-infiltration static constexpr int switch1Idx = Indices::switch1Idx; // wetting component switch index static constexpr int switch2Idx = Indices::switch2Idx; // Supposedly this is the index for the contaminant... static constexpr int temperatureIdx = Indices::temperatureIdx; // equation indices // inspired from Test problem infiltration in 3p3c static constexpr int contiWEqIdx = Indices::conti0EqIdx + FluidSystem::wCompIdx; //!< Index of the mass conservation equation for the water component static constexpr int contiNEqIdx = Indices::conti0EqIdx + FluidSystem::nCompIdx; //!< Index of the mass conservation equation for the contaminant component static constexpr int contiGEqIdx = Indices::conti0EqIdx + FluidSystem::gCompIdx; //!< Index of the mass conservation equation for the gas component static constexpr int energyEqIdx = Indices::energyEqIdx; static constexpr int H2OIdx = FluidSystem::H2OIdx; static constexpr int AirIdx = FluidSystem::AirIdx; static constexpr int NAPLIdx = FluidSystem::NAPLIdx; // phase indices static constexpr int wPhaseIdx = FluidSystem::wPhaseIdx; static constexpr int gPhaseIdx = FluidSystem::gPhaseIdx; static constexpr int nPhaseIdx = FluidSystem::nPhaseIdx; // phase states static constexpr int threePhases = Indices::threePhases; static constexpr int wgPhaseOnly = Indices::wgPhaseOnly; static constexpr int gnPhaseOnly = Indices::gnPhaseOnly; static constexpr int wnPhaseOnly = Indices::wnPhaseOnly; static constexpr int gPhaseOnly = Indices::gPhaseOnly; static constexpr int wPhaseOnly = Indices::wPhaseOnly; using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>; using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>; using Element = typename GridView::template Codim<0>::Entity; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView; using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>; using Grid = GetPropType<TypeTag, Properties::Grid>; using SubControlVolume = typename FVElementGeometry::SubControlVolume; using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView; using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView; using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; using EnergyLocalResidual = GetPropType<TypeTag, Properties::EnergyLocalResidual>; using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>; using MaterialLaw = typename HeatSourceProblem::SpatialParams::MaterialLaw; using MaterialLawParams = typename MaterialLaw::Params; using ThermalConductivityModel = GetPropType<TypeTag, Properties::ThermalConductivityModel>; static constexpr auto discMethod = GetPropType<TypeTag, Properties::FVGridGeometry>::discMethod; //! Property that defines whether mole or mass fractions are used static constexpr bool useMoles = ModelTraits::useMoles(); public: // this value will be used for comparisons static constexpr Scalar eps_ = 1e-4; // boundaries from gmsh IDs enum { diagId = 1, topId = 2, bottomId = 3, frontId = 4, leftId = 5, saeBottomId = 6, saeSideId = 7, hsId = 8, }; // // Also some types will be necessary. HeatSourceProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<const GridData<Grid>> gridData) : ParentType(fvGridGeometry), gridData_(gridData), heatSource_() { initialTemperature_ = getParam<Scalar>("Problem.InitialTemperature"); initialWaterSaturation_ = getParam<Scalar>("Problem.InitialWaterSaturation"); //initialNAPLSaturation_ = getParam<Scalar>("Problem.InitialNAPLSaturation"); saePressure_ = getParam<Scalar>("Problem.SaePressure"); saeRadius_ = getParam<Scalar>("Problem.SaeRadius"); saeBottom_ = getParam<Scalar>("Problem.SaeBottom"); saeHeight_ = getParam<Scalar>("Problem.SaeHeight"); saeBoundaryThickness_ = getParam<Scalar>("Problem.SaeBoundaryThickness"); // vtk output initialization effThermalConductivity_.resize(this->fvGridGeometry().numDofs()); liquidEnthalpy_.resize(this->fvGridGeometry().numDofs()); gasEnthalpy_.resize(this->fvGridGeometry().numDofs()); // set up solid and fluid system FluidSystem::init(); // print info about setup // this->info(); name_ = getParam<std::string>("Problem.Name"); // precompute the boundary types for the box method from the cell-centered boundary types scvfToScvBoundaryTypes_.computeBoundaryTypes(*this); } /****************************************************************************************** * Boundaries, Initial Values & Sources ******************************************************************************************/ /* Due to symmetry, only a wedge of the domain containing half a heat source * is simulated, and the hypotenuse side and the side adjacent to the * extraction well are symmetry planes, and therefore no-flow * boundaries. The extraction well and the heat source are cut out of the * domain and are implemented as solution dependent von Neumann * boundaries. The heat source boundaries is a no-flow boundary for air and * water, and the heat flux is calculated as if the heat source boundary is * at the given temperature, limited by the maximum heat flux (both * parameters can be set in the input-file). At the soil air extraction well * boundary, heat, water, and air component fluxes are calculated based on * the air phase flux, that is dependent on the pressure difference between * domain and extraction well (also given in the input file). Top, bottom, * and back side are no-flow boundaries. * * The boundaries are assigned using boundary markers from gmsh, with the * following meanings: * * 1 : diagonal side * 2 : top * 3 : bottom * 4 : front * 5 : left * 6 : soil air extraction bottom * 7 : soil air extraction mantle * 8 : heat source */ /*! * \fn * \brief The boundary condition types */ BoundaryTypes boundaryTypes(const Element &element, const SubControlVolume &scv) const { return scvfToScvBoundaryTypes_.boundaryTypes(scv); } BoundaryTypes boundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const { BoundaryTypes bcTypes; const auto boundaryMarkerId = gridData_->getBoundaryDomainMarker(scvf.boundaryFlag()); if (boundaryMarkerId == bottomId) bcTypes.setAllDirichlet(); else bcTypes.setAllNeumann(); return bcTypes; } /* * Solution dependent von Neumann boundaries * * We need special treatment of the boundaries of the soil air extraction * well and at the top and the bottom, since in these cases we want to have * heat and water transported with the air flow. * * At the sides, we just use normal no-flow */ NumEqVector neumann(const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, const ElementFluxVariablesCache& elemFluxVarsCache, const SubControlVolumeFace& scvf) const { NumEqVector flux(0.0); const auto& ipGlobal = scvf.ipGlobal(); // In order to find the boundary marker ID of the scvf for which this // method was called, we need to find the intersection this scvf belongs to const auto boundaryMarkerId = gridData_->getBoundaryDomainMarker(scvf.boundaryFlag()); // No flow boundaries // TODO: Top also to be handles as no flow boundary?! if (boundaryMarkerId == diagId || boundaryMarkerId == frontId || boundaryMarkerId == topId || boundaryMarkerId == leftId) { return flux; } else if (boundaryMarkerId == hsId) { // Heat source boundary: only energy input // // Implemented vapor injection here: disactivated for simplicity reasons //flux[contiGEqIdx] = -0.2; //flux[contiWEqIdx] = -0.1435; //flux[contiNEqIdx] = 0.0; // kein NAPL einpumpen // //if (!heatSource_.isHeated(ipGlobal[2])) return flux; // get heat flux from temperature inside the domain auto sol = elementSolution(element, elemVolVars, fvGeometry); const Scalar T = sol[fvGeometry.scv(scvf.insideScvIdx()).localDofIndex()][temperatureIdx]; flux[contiGEqIdx] = 0.0; flux[contiWEqIdx] = 0.0; flux[contiNEqIdx] = 0.0; flux[energyEqIdx] = heatSource_.heatFlux(T); return flux; } else if (boundaryMarkerId == saeBottomId || boundaryMarkerId == saeSideId) { return flux; // TODO: remove this when program works properly // upper part is no-flow if (ipGlobal[2] > saeBottom_ + saeHeight_) return flux; // get volume variables const auto& volVars = elemVolVars[scvf.insideScvIdx()]; // get pressure gradient const Scalar pa = volVars.pressure(gPhaseIdx); Scalar gasPhaseDarcyFlux = (pa - saePressure_) / saeBoundaryThickness_; // if the pressure in the cell is higher than the one set on the // boundary, the flux should be outwards, i.e. positive. // since unitOuterNormal and pressureGradient here both point // outwards, there should be no negative sign in front of the // permeability // assuming relative permeability of wall = 1 const Scalar mu = volVars.fluidState().viscosity(gPhaseIdx); gasPhaseDarcyFlux *= volVars.permeability() / mu; // the component flux over the boundary is what is carried away in the air flux const Scalar gasMolarDensity = volVars.molarDensity(gPhaseIdx); const Scalar airFraction = volVars.moleFraction(gPhaseIdx, AirIdx); const Scalar H2OFraction = volVars.moleFraction(gPhaseIdx, H2OIdx); const Scalar NAPLFraction = volVars.moleFraction(gPhaseIdx, NAPLIdx); // transport of air and water in gas phase flux[contiGEqIdx] = gasMolarDensity * airFraction * gasPhaseDarcyFlux; flux[contiWEqIdx] = gasMolarDensity * H2OFraction * gasPhaseDarcyFlux; flux[contiNEqIdx] = gasMolarDensity * NAPLFraction * gasPhaseDarcyFlux; // transport of latent heat const Scalar gasEnthalpy = volVars.enthalpy(gPhaseIdx); //[J/kg] const Scalar gasMassDensity = volVars.density(gPhaseIdx); flux[energyEqIdx] = gasMassDensity * gasEnthalpy * gasPhaseDarcyFlux; return flux; } else { std::cerr << "no such boundary marker id!" << std::endl; } return flux; } NumEqVector source(const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, const SubControlVolume& scv) const { NumEqVector values(0.0); return values; } PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const { // this is only called for the boundaries that are similar to the // initial condition return initialAtPos(globalPos); } PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const { PrimaryVariables priVars (0.0); priVars[switch1Idx] = initialWaterSaturation_; priVars[switch2Idx] = 0;// initialNAPLSaturation_; priVars[temperatureIdx] = initialTemperature_; const Scalar pa0 = 101325; //wikipedia //const Scalar densityAir = 1.2466; //from Wikipedia const Scalar densityWater = 1e3; // very rough, just for the beginning! const Scalar g = 9.81; // const Scalar pa = pa0 - globalPos[dimWorld-1] * densityAir * g; // disactivated for simplicity reasons priVars[pressureIdx] = pa0; // Over the water table we assume everywhere atmosphere pressure priVars.setState(wgPhaseOnly); // Vadose zone // Enplacement of NAPL contamination; disactivated for simplicity reasons /* if (isInContaminationZone(globalPos)) { priVars.setState(threePhases); priVars[switch2Idx] = 0.07; } // else // priVars.setState(wgPhaseOnly); */ // from 0 to 2.5 height I will place a fully saturated zone, upwards it is a vadose zone with p=p_atm /* // disactivated for simplicity reasons if( globalPos[2]< 2.5 +eps_) { priVars.setState(wPhaseOnly); // Here is only water. // TODO: Take water density from material law! priVars[pressureIdx] = pa0 + g*densityWater *(2.5- globalPos[2]); //priVars[switch1Idx] = 0.98; } */ // TODO: The boundary pressure is to be placed here. return priVars; } /**************************************************************************** * Function to be used in the time loop ****************************************************************************/ /****************************************************************************************** * Getters ******************************************************************************************/ const std::string& name() const { return name_; } const Scalar saeBottom() const { return saeBottom_; } const Scalar saeTop() const { return saeBottom_ + saeHeight_; } const Scalar saeRadius() const { return saeRadius_; } /****************************************************************************************** * Plots ******************************************************************************************/ void plotWaterDensity() const { GnuplotInterface<Scalar> gnuplot(true); gnuplot.setOpenPlotWindow(true); std::vector<Scalar> T; std::vector<Scalar> rho_l, rho_g; for (double temp = 283.15; temp < 283.15 + 700; temp += 1.0) { T.push_back(temp - 273.15); rho_l.push_back(Components::H2O<Scalar>::liquidDensity(temp, 1e5)); rho_g.push_back(Components::H2O<Scalar>::gasDensity(temp, 1e5)); } gnuplot.setXlabel("Temperature [°C]"); gnuplot.setYlabel("Liquid Density of Water [kg/m**3]"); gnuplot.addDataSetToPlot(T, rho_l, "T-rho_l", "w l"); gnuplot.setOption("set grid"); gnuplot.plot("T-rho_l"); gnuplot.resetAll(); gnuplot.setOpenPlotWindow(true); gnuplot.setXlabel("Temperature [°C]"); gnuplot.setYlabel("Gas Density of Water [kg/m**3]"); gnuplot.addDataSetToPlot(T, rho_g, "T-rho_g", "w l"); gnuplot.setOption("set grid"); gnuplot.plot("T-rho_g"); } void plotWaterRetentionCurve() const { GnuplotInterface<Scalar> gnuplot(true); gnuplot.resetAll(); gnuplot.setOpenPlotWindow(true); std::vector<Scalar> Sw; std::vector<Scalar> pc; for (Scalar sw = 0.0; sw <= 1.0; sw += 0.01) { Sw.push_back(sw); pc.push_back(MaterialLaw::pc(this->spatialParams().materialLawParamsAtPos({0.0, 0.0, 0.0}), sw)); } gnuplot.setXlabel("Pressure"); gnuplot.setYlabel("Water Saturation"); gnuplot.addDataSetToPlot(Sw, pc, "Sw-pc", "w l"); gnuplot.setOption("set grid"); gnuplot.plot("Sw-pc"); } private: // checks whether a point is located inside the contamination zone bool isInContaminationZone(const GlobalPosition &globalPos) const { // TODO: Does it work well? Still not able to see // Place a small NAPL-contamination at : // -3 <= x <= -2 ; 2<=y <=3 ; 4<=z<= 5 return ( Dune::FloatCmp::ge<Scalar>(globalPos[0], -3) && Dune::FloatCmp::le<Scalar>(globalPos[0], -2) && Dune::FloatCmp::ge<Scalar>(globalPos[1], 2 ) && Dune::FloatCmp::le<Scalar>(globalPos[1], 3 ) && Dune::FloatCmp::ge<Scalar>(globalPos[2], 4 ) && Dune::FloatCmp::le<Scalar>(globalPos[2], 5 ) ); } /****************************************************************************************** * Initialization ******************************************************************************************/ /****************************************************************************************** * Domain definitions ******************************************************************************************/ // bool isFront(const GlobalPosition &pos) const { return pos[0] < eps_; } // bool isBack(const GlobalPosition &pos) const { return pos[0] > domainSize_[0] - eps_; } // bool isLeft(const GlobalPosition &pos) const { return pos[1] < eps_; } // bool isRight(const GlobalPosition &pos) const { return pos[1] > domainSize_[1] - eps_; } // bool isBottom(const GlobalPosition &pos) const { return pos[2] < eps_; } // bool isTop(const GlobalPosition &pos) const { return pos[2] > domainSize_[2] - eps_; } // bool isHe(const GlobalPosition& pos) const // { // return pow(hePosition_[0] - pos[0], 2) + pow(hePosition_[1] - pos[1], 2) <= pow(1.1*heRadius_, 2); // } // bool isSae(const GlobalPosition& pos) const // { // return pow(pos[0], 2) + pow(pos[1], 2) <= pow(1.1*saeRadius_, 2); // } /****************************************************************************************** * Member variables *******************************************************************************************/ std::shared_ptr<const GridData<Grid>> gridData_; HeatSource<TypeTag> heatSource_; Scalar initialTemperature_; Scalar initialWaterSaturation_, initialNAPLSaturation_; Scalar saeRadius_; Scalar saePressure_; Scalar saeHeight_; Scalar saeBottom_; Scalar saeBoundaryThickness_; GlobalPosition domainSize_; // output bool plotFluidMatrixInteractions_; std::vector<Scalar> effThermalConductivity_; std::vector<Scalar> liquidEnthalpy_; std::vector<Scalar> gasEnthalpy_; std::string name_; ScvfToScvBoundaryTypes<BoundaryTypes, discMethod> scvfToScvBoundaryTypes_; }; } // 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/>. * *****************************************************************************/ #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 <dune/istl/io.hh> #include <dumux/discretization/method.hh> #include <dumux/common/properties.hh> #include <dumux/common/parameters.hh> #include <dumux/common/dumuxmessage.hh> #include <dumux/common/defaultusagemessage.hh> #include <dumux/common/timeloop.hh> #include <dumux/linear/amgbackend.hh> #include <dumux/linear/seqsolverbackend.hh> #include <dumux/nonlinear/newtonsolver.hh> #include <dumux/assembly/fvassembler.hh> #include <dumux/assembly/diffmethod.hh> #include <dumux/io/vtkoutputmodule.hh> #include <dumux/io/grid/gridmanager.hh> #include <dumux/nonlinear/newtonconvergencewriter.hh> #include "trackedquantities.hh" #include "tracker.hh" #include "problem.hh" int main(int argc, char** argv) try { using namespace Dumux; // define the type tag for this problem using TypeTag = Properties::TTag::HeatSourceBox; //////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////// // 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); // initialize parameter tree Parameters::init(argc, argv, "params.input"); ////////////////////////////////////////////////////////////////////// // try to create a grid (from the given grid file or the input file) ///////////////////////////////////////////////////////////////////// using Grid = GetPropType<TypeTag, Properties::Grid>; GridManager<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::FVGridGeometry>; 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, gridManager.getGridData()); // the solution vector using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; SolutionVector x(fvGridGeometry->numDofs()); problem->applyInitialSolution(x); auto xOld = x; // the grid variables using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry); gridVariables->init(x); // intialize the vtk output module using IOFields = GetPropType<TypeTag, Properties::IOFields>; VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name()); using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>; vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables)); IOFields::initOutputModule(vtkWriter); //!< Add model specific output fields // vtkWriter.addField(problem->getEffThermalConductivity(), "lambda_eff"); // vtkWriter.addField(problem->getLiquidEnthalpy(), "h_liquid"); // vtkWriter.addField(problem->getGasEnthalpy(), "h_gas"); // problem->updateVtkOutput(x); vtkWriter.write(0.0); // get some time loop parameters using Scalar = GetPropType<TypeTag, Properties::Scalar>; const Scalar tEnd = getParam<Scalar>("TimeLoop.TEnd"); Scalar dt = getParam<Scalar>("TimeLoop.DtInitial"); const Scalar dtMax = getParam<Scalar>("TimeLoop.DtMax"); // instantiate time loop auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd); timeLoop->setMaxTimeStepSize(dtMax); // set some check points for the time loop const Scalar dtOut = getParam<Scalar>("TimeLoop.DtOut"); timeLoop->setPeriodicCheckPoint(dtOut); // the assembler with time loop for instationary problem using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>; auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop); // the linear solver using LinearSolver = AMGBackend<TypeTag>; auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper()); //using LinearSolver = ILU0BiCGSTABBackend; //auto linearSolver = std::make_shared<LinearSolver>(); // the non-linear solver using NewtonSolver = NewtonSolver<Assembler, LinearSolver>; NewtonSolver nonLinearSolver(assembler, linearSolver); //the convergence writer //bool writeConvergence = getParam<bool>("Newton.WriteConvergence"); //using GridView = GetPropType<TypeTag, Properties::GridView>; //using NewtonConvergenceWriter = Dumux::NewtonConvergenceWriter<GridView, SolutionVector>; //auto convergenceWriter = std::make_shared<NewtonConvergenceWriter>(leafGridView, fvGridGeometry->numDofs()); //the convergence writer using NewtonConvergenceWriter = Dumux::NewtonConvergenceWriter<FVGridGeometry, SolutionVector>; auto convergenceWriter = std::make_shared<NewtonConvergenceWriter>(*fvGridGeometry); nonLinearSolver.attachConvergenceWriter(convergenceWriter); // problem->plotWaterDensity(); // problem->plotWaterRetentionCurve(); // // Tracking of different quantities // using GridView = GetPropType<TypeTag, Properties::GridView>; // using Element = typename GridView::template Codim<0>::Entity; // using GlobalPosition = typename Element::Geometry::GlobalCoordinate; // using namespace DumuxVegasGeothermal::Tracking; // // auto trackedFluxQuantities = std::make_tuple( // FluxAcrossMarkedInterface<TypeTag>(problem->saeSideId, gridManager.getGridData()), // FluxAcrossMarkedInterface<TypeTag>(problem->saeBottomId, gridManager.getGridData()) // ); // auto trackedVolumeQuantities = std::make_tuple( // // track temperature at soil air extraction well // TrackedTemperatureAtPos<TypeTag>( // "T_sae_bot", GlobalPosition({-problem->saeRadius() - problem->eps_, // problem->eps_, // problem->saeBottom()})), // TrackedTemperatureAtPos<TypeTag>( // "T_sae_top", GlobalPosition({-problem->saeRadius() - problem->eps_, // problem->eps_, // problem->saeTop()})) // ); // using TrackedFluxQuantities = decltype(trackedFluxQuantities); // using TrackedVolumeQuantities = decltype(trackedVolumeQuantities); // using Tracker = Tracker<TypeTag, TrackedFluxQuantities, TrackedVolumeQuantities>; // auto tracker = std::make_shared<Tracker>(trackedFluxQuantities, trackedVolumeQuantities); // tracker->initializeOutputFile(); // time loop timeLoop->start(); do { // set previous solution for storage evaluations assembler->setPreviousSolution(xOld); // solve the non-linear system with time step control nonLinearSolver.solve(x, *timeLoop); // make the new solution the old solution xOld = x; gridVariables->advanceTimeStep(); // based on how easy it was to solve for the current // timestep, we can suggest a new timestep: dt = nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()); // advance to the time loop to the next step // In case a checkpoint is ahead (i.e. in the next loop) this already // sets the time step accordingly timeLoop->advanceTimeStep(); // // write out tracked variables // tracker->updateTrackedQuantities(*problem, *gridVariables, x, *assembler); // tracker->appendToOutputFile(timeLoop->time()); // write vtk output at checkpoint if (timeLoop->isCheckPoint()) { // problem->updateVtkOutput(x); vtkWriter.write(timeLoop->time()); } // report statistics of this time step timeLoop->reportTimeStep(); timeLoop->setTimeStepSize(dt); } while (!timeLoop->finished()); timeLoop->finalize(leafGridView.comm()); //////////////////////////////////////////////////////////// // finalize, print dumux message to say goodbye //////////////////////////////////////////////////////////// // print dumux end message if (mpiHelper.rank() == 0) DumuxMessage::print(/*firstCall=*/false); return 0; } catch (Dumux::ParameterException &e) { std::cerr << std::endl << e << " ---> Abort!" << std::endl; return 1; } catch (Dune::DGFException & e) { std::cerr << "DGF exception thrown (" << e << "). Most likely, the DGF file name is wrong " "or the DGF file is corrupted, " "e.g. missing hash at end of file or wrong number (dimensions) of entries." << " ---> Abort!" << std::endl; return 2; } catch (Dune::Exception &e) { std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl; return 3; } catch (...) { std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; return 4; }
params.input
Description: Binary data
// -*- 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 3 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 * \ingroup VegasGeothermal * \brief Definition of the spatial parameters for the Vegas container. */ #ifndef DUMUX_VEGAS_GEOTHERMAL_REALWORLD_SPATIAL_PARAMS_HH #define DUMUX_VEGAS_GEOTHERMAL_REALWORLD_SPATIAL_PARAMS_HH #include <iostream> #include <dumux/io/gnuplotinterface.hh> #include <dumux/io/ploteffectivediffusivitymodel.hh> #include <dumux/io/plotmateriallaw.hh> #include <dumux/io/plotthermalconductivitymodel.hh> #include <dumux/porousmediumflow/properties.hh> #include <dumux/material/spatialparams/fv.hh> #include <dumux/material/fluidmatrixinteractions/3p/regularizedparkervangen3p.hh> //#include <dumux/material/fluidmatrixinteractions/3p/vangenuchtenoftemperature.hh> // TODO: Write this for 3p! Only if necessary! // From remscen: #include <dune/common/float_cmp.hh> #include <dumux/porousmediumflow/properties.hh> #include <dumux/material/spatialparams/fv.hh> #include <dumux/material/fluidmatrixinteractions/3p/regularizedparkervangen3p.hh> #include <dumux/material/fluidmatrixinteractions/3p/regularizedparkervangen3pparams.hh> #include <dumux/material/fluidmatrixinteractions/3p/efftoabslaw.hh> #include <dumux/material/fluidmatrixinteractions/3p/efftoabslaw.hh> #include <dumux/material/solidstates/inertsolidstate.hh> #include <dumux/material/solidsystems/inertsolidphase.hh> namespace DumuxVegasGeothermal { using namespace Dumux; /*! * \ingroup VegasGeothermal * \brief Definition of the spatial parameters for the Vegas container. */ template<class FVGridGeometry, class Scalar> class HorizontallyLayeredSpatialParams : public FVSpatialParams<FVGridGeometry, Scalar, HorizontallyLayeredSpatialParams<FVGridGeometry, Scalar>> { using Grid = typename FVGridGeometry::Grid; 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 = FVSpatialParams<FVGridGeometry, Scalar, HorizontallyLayeredSpatialParams<FVGridGeometry, Scalar>>; static constexpr int dimWorld = GridView::dimensionworld; using EffectiveLaw = ParkerVanGen3P<Scalar>; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; public: // this value will be used for comparisons static constexpr Scalar eps_ = 1e-4; //! Export the type used for the permeability using PermeabilityType = Scalar; //! Export the type used for the material law using MaterialLaw = EffToAbsLaw<EffectiveLaw>; using MaterialLawParams = typename MaterialLaw::Params; HorizontallyLayeredSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry) : ParentType(fvGridGeometry) { // read soil layering soiltop_ = getParam<std::vector<Scalar>>("SoilLayering.ZTop"); porosity_ = getParam<std::vector<Scalar>>("SoilLayering.Porosity"); permeability_ = getParam<std::vector<Scalar>>("SoilLayering.Permeability"); solidDensity_ = getParam<std::vector<Scalar>>("SoilLayering.SolidDensity"); solidHeatCapacity_ = getParam<std::vector<Scalar>>("SoilLayering.SolidHeatCapacity"); solidThermalConductivity_ = getParam<std::vector<Scalar>>("SoilLayering.SolidThermalConductivity"); swr_ = getParam<std::vector<Scalar>>("SoilLayering.Swr"); snr_ = getParam<std::vector<Scalar>>("SoilLayering.Snr"); vgAlpha_ = getParam<std::vector<Scalar>>("SoilLayering.VgAlpha"); vgn_ = getParam<std::vector<Scalar>>("SoilLayering.Vgn"); pcLowSw_ = getParam<std::vector<Scalar>>("SoilLayering.PcLowSw"); pcHighSw_ = getParam<std::vector<Scalar>>("SoilLayering.PcHighSw"); numLayers_ = soiltop_.size(); // material law params materialLawParams_.resize(numLayers_); // intrinsic permeabilites (taken from remediationscenarios fineK_ = 5.47e-10; coarseK_ = 9.14e-10; for (std::size_t i = 0; i < numLayers_; i++) { materialLawParams_[i].setSwr(swr_[i]); materialLawParams_[i].setSnr(snr_[i]); materialLawParams_[i].setVgAlpha(vgAlpha_[i]); materialLawParams_[i].setVgn(vgn_[i]); } } // Place different permeabilites Scalar permeabilityAtPos(const GlobalPosition& globalPos) const { if(globalPos[2] > 7.5 -eps_) { // return fineK_[getSoilIdx_(globalPos)]; return fineK_; } return coarseK_; //return coarseK_[getSoilIdx_(globalPos)]; } Scalar porosityAtPos(const GlobalPosition& globalPos) const { return porosity_[getSoilIdx_(globalPos)]; } const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const { return materialLawParams_[getSoilIdx_(globalPos)]; } // For hydrostatic pressure see dumux-course-problem slide 12 template<class FluidSystem> int wettingPhaseAtPos(const GlobalPosition& globalPos) const { return FluidSystem::H2OIdx; } template<class SolidSystem> Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos, int compIdx) const { return 0; } protected: std::size_t getSoilIdx_(const GlobalPosition& pos) const { const auto z = pos[dimWorld-1]; for (std::size_t i = 0; i < numLayers_; i++) { if (z <= soiltop_[i]) return i; } DUNE_THROW(Dune::InvalidStateException, "No soil definition at z = " + std::to_string(z)); } Scalar fineK_; Scalar coarseK_; std::size_t numLayers_; std::vector<Scalar> soiltop_; std::vector<Scalar> porosity_; std::vector<Scalar> permeability_; std::vector<Scalar> solidDensity_; std::vector<Scalar> solidHeatCapacity_; std::vector<Scalar> solidThermalConductivity_; std::vector<Scalar> swr_; std::vector<Scalar> snr_; std::vector<Scalar> vgAlpha_; std::vector<Scalar> vgn_; std::vector<Scalar> pcLowSw_; std::vector<Scalar> pcHighSw_; std::vector<MaterialLawParams> materialLawParams_; }; } // end namespace Dumux #endif
_______________________________________________ Dumux mailing list Dumux@listserv.uni-stuttgart.de https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
