diff --git a/test/mixeddimension/embedded/1p_richards/richardstestproblem.hh b/test/mixeddimension/embedded/1p_richards/richardstestproblem.hh index de02aaad8565ffcc078905e35c73213b9ceb6bb9..d8a1b4d5b796ea6daaab9c28cd569c4c9c03ef7c 100644 --- a/test/mixeddimension/embedded/1p_richards/richardstestproblem.hh +++ b/test/mixeddimension/embedded/1p_richards/richardstestproblem.hh @@ -232,6 +232,37 @@ public: source = sourceValue*source.quadratureWeight()*source.integrationElement(); } + //! Called after every time step + //! Output the total global exchange term + void postTimeStep() + { + ParentType::postTimeStep(); + + PrimaryVariables source(0.0); + + if (!(this->timeManager().time() < 0.0)) + { + for (const auto& element : elements(this->gridView())) + { + auto fvGeometry = localView(this->model().globalFvGeometry()); + fvGeometry.bindElement(element); + + auto elemVolVars = localView(this->model().curGlobalVolVars()); + elemVolVars.bindElement(element, fvGeometry, this->model().curSol()); + + for (auto&& scv : scvs(fvGeometry)) + { + auto pointSources = this->scvPointSources(element, fvGeometry, elemVolVars, scv); + pointSources *= scv.volume()*elemVolVars[scv].extrusionFactor(); + source += pointSources; + } + } + } + + std::cout << "Global integrated source (soil): " << source << " (kg/s) / " + << source*3600*24*1000 << " (g/day)" << '\n'; + } + // \} /*! * \name Boundary conditions diff --git a/test/mixeddimension/embedded/1p_richards/rootsystemtestproblem.hh b/test/mixeddimension/embedded/1p_richards/rootsystemtestproblem.hh index 6507c30fa3deb7665dd150f24cf5d1cec4bd1150..5ed1fb529330a351b1ba8da9f46b20688fe5c324 100644 --- a/test/mixeddimension/embedded/1p_richards/rootsystemtestproblem.hh +++ b/test/mixeddimension/embedded/1p_richards/rootsystemtestproblem.hh @@ -59,6 +59,13 @@ public: // Set the grid type SET_TYPE_PROP(RootsystemTestProblem, Grid, Dune::FoamGrid<1, 3>); +SET_BOOL_PROP(RootsystemTestProblem, EnableGlobalFVGeometryCache, true); +SET_BOOL_PROP(RootsystemTestProblem, EnableGlobalVolumeVariablesCache, true); +SET_BOOL_PROP(RootsystemTestProblem, EnableGlobalFluxVariablesCache, true); +SET_BOOL_PROP(RootsystemTestProblem, SolutionDependentAdvection, false); +SET_BOOL_PROP(RootsystemTestProblem, SolutionDependentMolecularDiffusion, false); +SET_BOOL_PROP(RootsystemTestProblem, SolutionDependentHeatConduction, false); + // Set the problem property SET_TYPE_PROP(RootsystemTestProblem, Problem, RootsystemTestProblem<TypeTag>); @@ -213,12 +220,12 @@ public: const SubControlVolumeFace& scvf) const { PrimaryVariables values(0.0); - if (scvf.center()[2] + eps_ > this->bBoxMax()[2] ) + if (scvf.center()[2] + eps_ > this->bBoxMax()[2]) { const auto r = this->spatialParams().radius(this->gridView().indexSet().index(element)); values[conti0EqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, - BoundaryConditions.TranspirationRate)/(M_PI*r*r); + BoundaryConditions.TranspirationRate)/(M_PI*r*r)/scvf.area(); } return values; @@ -298,6 +305,37 @@ public: source = sourceValue*source.quadratureWeight()*source.integrationElement(); } + //! Called after every time step + //! Output the total global exchange term + void postTimeStep() + { + ParentType::postTimeStep(); + + PrimaryVariables source(0.0); + + if (!(this->timeManager().time() < 0.0)) + { + for (const auto& element : elements(this->gridView())) + { + auto fvGeometry = localView(this->model().globalFvGeometry()); + fvGeometry.bindElement(element); + + auto elemVolVars = localView(this->model().curGlobalVolVars()); + elemVolVars.bindElement(element, fvGeometry, this->model().curSol()); + + for (auto&& scv : scvs(fvGeometry)) + { + auto pointSources = this->scvPointSources(element, fvGeometry, elemVolVars, scv); + pointSources *= scv.volume()*elemVolVars[scv].extrusionFactor(); + source += pointSources; + } + } + } + + std::cout << "Global integrated source (root): " << source << " (kg/s) / " + << source*3600*24*1000 << " (g/day)" << '\n'; + } + bool shouldWriteRestartFile() const { return false; diff --git a/test/mixeddimension/embedded/1p_richards/rootsystemtestspatialparams.hh b/test/mixeddimension/embedded/1p_richards/rootsystemtestspatialparams.hh index 2b1873181b139af6d7383895d4291506cac860dc..1bb9adf3876ef844d8edbd926084e1261420b205 100644 --- a/test/mixeddimension/embedded/1p_richards/rootsystemtestspatialparams.hh +++ b/test/mixeddimension/embedded/1p_richards/rootsystemtestspatialparams.hh @@ -25,6 +25,7 @@ #define DUMUX_ROOTSYSTEM_TEST_SPATIALPARAMS_HH #include <dumux/material/spatialparams/implicit1p.hh> +#include <dumux/material/components/simpleh2o.hh> namespace Dumux { @@ -45,7 +46,9 @@ class RootsystemTestSpatialParams: public ImplicitSpatialParamsOneP<TypeTag> using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using Element = typename GridView::template Codim<0>::Entity; using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); + using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector); + using Indices = typename GET_PROP_TYPE(TypeTag, Indices); enum { // Grid and world dimension dim = GridView::dimension, @@ -102,8 +105,10 @@ public: rootParams_[eIdx].radialPerm = Kr_; //Kr } else //order >= 3 + { rootParams_[eIdx].axialPerm = Kx_; //Kx rootParams_[eIdx].radialPerm = Kr_; //Kr + } } } @@ -111,15 +116,15 @@ public: * \brief Return the intrinsic permeability for the current sub-control volume in [m^2]. * * \param ipGlobal The integration point - * \note Kx has units [m^4] so we have to divide by the cross-section area + * \note Kx has units [m^4/(Pa*s)] so we have to divide by the cross-section area + * and multiply with a characteristic viscosity */ - Scalar permeability(const Element& element, - const SubControlVolume& scv, - const ElementSolutionVector& elemSol) const + PermeabilityType permeability(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const { - auto eIdx = gridView_.indexSet().index(element); - const auto r = radius(eIdx); - return Kx_ / (M_PI*r*r); + const Scalar r = rootParams(element).radius; + return Kx_ / (M_PI*r*r) * SimpleH2O<Scalar>::liquidViscosity(285.15, elemSol[0][Indices::pressureIdx]); } /*! diff --git a/test/mixeddimension/embedded/1p_richards/rositestproblem.hh b/test/mixeddimension/embedded/1p_richards/rositestproblem.hh index ef5df270c8ec79da531f8dcc8330afb4c76225d7..4789371b44fd035d54aac44efa251f9fd2e836c9 100644 --- a/test/mixeddimension/embedded/1p_richards/rositestproblem.hh +++ b/test/mixeddimension/embedded/1p_richards/rositestproblem.hh @@ -29,6 +29,7 @@ #include <dumux/mixeddimension/problem.hh> #include <dumux/mixeddimension/embedded/cellcentered/bboxtreecouplingmanager.hh> +#include <dumux/mixeddimension/embedded/cellcentered/bboxtreecouplingmanagersimple.hh> #include <dumux/mixeddimension/integrationpointsource.hh> namespace Dumux @@ -44,7 +45,8 @@ NEW_TYPE_TAG(RosiTestProblem, INHERITS_FROM(MixedDimension)); SET_TYPE_PROP(RosiTestProblem, Problem, Dumux::RosiTestProblem<TypeTag>); // Set the coupling manager -SET_TYPE_PROP(RosiTestProblem, CouplingManager, Dumux::CCBBoxTreeEmbeddedCouplingManager<TypeTag>); +//SET_TYPE_PROP(RosiTestProblem, CouplingManager, Dumux::CCBBoxTreeEmbeddedCouplingManager<TypeTag>); +SET_TYPE_PROP(RosiTestProblem, CouplingManager, Dumux::CCBBoxTreeEmbeddedCouplingManagerSimple<TypeTag>); // Set the two sub-problems of the global problem SET_TYPE_PROP(RosiTestProblem, LowDimProblemTypeTag, TTAG(RootsystemTestCCProblem)); diff --git a/test/mixeddimension/embedded/1p_richards/test_rosi.input b/test/mixeddimension/embedded/1p_richards/test_rosi.input index c0dbe86c18e2c7497af1591fc1907621aa2fe264..4835a18d195c2e0489138defcce7ab6a2e331c7c 100644 --- a/test/mixeddimension/embedded/1p_richards/test_rosi.input +++ b/test/mixeddimension/embedded/1p_richards/test_rosi.input @@ -18,7 +18,6 @@ Cells = 20 20 20 [Problem] Name = rosi -EnableGravity = false [SpatialParams] Permeability = 2.57e-12 # [m^2]