| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /** | ||
| 8 | * \file | ||
| 9 | * \ingroup EmbeddedTests | ||
| 10 | * \brief Definition of a problem, for the 1p2c problem: | ||
| 11 | * Component transport of oxygen in interstitial fluid. | ||
| 12 | */ | ||
| 13 | |||
| 14 | #ifndef DUMUX_TISSUE_PROBLEM_HH | ||
| 15 | #define DUMUX_TISSUE_PROBLEM_HH | ||
| 16 | |||
| 17 | #include <dune/geometry/quadraturerules.hh> | ||
| 18 | #include <dune/localfunctions/lagrange/pqkfactory.hh> | ||
| 19 | |||
| 20 | #include <dumux/common/boundarytypes.hh> | ||
| 21 | #include <dumux/common/math.hh> | ||
| 22 | #include <dumux/common/parameters.hh> | ||
| 23 | #include <dumux/common/properties.hh> | ||
| 24 | #include <dumux/common/numeqvector.hh> | ||
| 25 | |||
| 26 | #include <dumux/porousmediumflow/problem.hh> | ||
| 27 | |||
| 28 | #include <dumux/multidomain/embedded/couplingmanager1d3d.hh> // for coupling mode | ||
| 29 | |||
| 30 | namespace Dumux { | ||
| 31 | |||
| 32 | /*! | ||
| 33 | * \ingroup EmbeddedTests | ||
| 34 | * \brief Definition of a problem, for the 1p2c problem: | ||
| 35 | * Component transport of oxygen in interstitial fluid. | ||
| 36 | */ | ||
| 37 | template <class TypeTag> | ||
| 38 | class TissueProblem : public PorousMediumFlowProblem<TypeTag> | ||
| 39 | { | ||
| 40 | using ParentType = PorousMediumFlowProblem<TypeTag>; | ||
| 41 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 42 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 43 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 44 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 45 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 46 | using GridView = typename GridGeometry::GridView; | ||
| 47 | using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; | ||
| 48 | using NumEqVector = Dumux::NumEqVector<PrimaryVariables>; | ||
| 49 | using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; | ||
| 50 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 51 | using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; | ||
| 52 | using PointSource = GetPropType<TypeTag, Properties::PointSource>; | ||
| 53 | using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; | ||
| 54 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 55 | using GlobalPosition = typename GridGeometry::GlobalCoordinate; | ||
| 56 | |||
| 57 | using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>; | ||
| 58 | |||
| 59 | public: | ||
| 60 | 19 | TissueProblem(std::shared_ptr<const GridGeometry> gridGeometry, | |
| 61 | std::shared_ptr<CouplingManager> couplingManager) | ||
| 62 | : ParentType(gridGeometry, "Tissue") | ||
| 63 |
5/18✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 19 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 19 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 19 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
|
57 | , couplingManager_(couplingManager) |
| 64 | { | ||
| 65 | // read parameters from input file | ||
| 66 |
9/26✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 19 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 19 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 19 times.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 19 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 19 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 19 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 19 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
|
38 | name_ = getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name"); |
| 67 | |||
| 68 |
3/6✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 19 times.
✗ Branch 8 not taken.
|
40 | exactPressure_.resize(this->gridGeometry().numDofs()); |
| 69 |
2/4✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19 times.
✗ Branch 5 not taken.
|
38 | auto fvGeometry = localView(this->gridGeometry()); |
| 70 |
4/8✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 19 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 931237 times.
✗ Branch 10 not taken.
|
1862493 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 71 | { | ||
| 72 | 931218 | fvGeometry.bindElement(element); | |
| 73 | |||
| 74 |
7/8✓ Branch 0 taken 109744 times.
✓ Branch 1 taken 13718 times.
✓ Branch 2 taken 1027244 times.
✓ Branch 3 taken 931218 times.
✓ Branch 4 taken 917500 times.
✓ Branch 5 taken 917500 times.
✓ Branch 7 taken 917500 times.
✗ Branch 8 not taken.
|
3916924 | for (auto&& scv : scvs(fvGeometry)) |
| 75 |
2/4✓ Branch 1 taken 1027244 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1027244 times.
✗ Branch 5 not taken.
|
2054488 | exactPressure_[scv.dofIndex()] = exactSolution(scv.dofPosition()); |
| 76 | } | ||
| 77 | 19 | } | |
| 78 | |||
| 79 | /*! | ||
| 80 | * \brief The problem name. | ||
| 81 | * | ||
| 82 | * This is used as a prefix for files generated by the simulation. | ||
| 83 | */ | ||
| 84 | const std::string& name() const | ||
| 85 |
1/2✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
|
19 | { return name_; } |
| 86 | |||
| 87 | // \} | ||
| 88 | |||
| 89 | /*! | ||
| 90 | * \name Boundary conditions | ||
| 91 | */ | ||
| 92 | // \{ | ||
| 93 | |||
| 94 | /*! | ||
| 95 | * \brief Specifies which kind of boundary condition should be | ||
| 96 | * used for which equation on a given boundary segment. | ||
| 97 | * | ||
| 98 | * \param globalPos The position for which the bc type should be evaluated | ||
| 99 | */ | ||
| 100 | BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const | ||
| 101 | { | ||
| 102 |
4/6✓ Branch 0 taken 186544 times.
✓ Branch 1 taken 57888 times.
✓ Branch 2 taken 173000 times.
✓ Branch 3 taken 55000 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
472432 | BoundaryTypes values; |
| 103 |
40/60✓ Branch 0 taken 186544 times.
✓ Branch 1 taken 57888 times.
✓ Branch 2 taken 186544 times.
✓ Branch 3 taken 57888 times.
✓ Branch 4 taken 186544 times.
✓ Branch 5 taken 57888 times.
✓ Branch 6 taken 186544 times.
✓ Branch 7 taken 57888 times.
✓ Branch 8 taken 186544 times.
✓ Branch 9 taken 57888 times.
✓ Branch 10 taken 57888 times.
✓ Branch 11 taken 128656 times.
✓ Branch 12 taken 57888 times.
✓ Branch 13 taken 128656 times.
✓ Branch 14 taken 57888 times.
✓ Branch 15 taken 128656 times.
✓ Branch 16 taken 57888 times.
✓ Branch 17 taken 128656 times.
✓ Branch 18 taken 57888 times.
✓ Branch 19 taken 128656 times.
✓ Branch 20 taken 173000 times.
✓ Branch 21 taken 55000 times.
✓ Branch 22 taken 173000 times.
✓ Branch 23 taken 55000 times.
✓ Branch 24 taken 173000 times.
✓ Branch 25 taken 55000 times.
✓ Branch 26 taken 173000 times.
✓ Branch 27 taken 55000 times.
✓ Branch 28 taken 173000 times.
✓ Branch 29 taken 55000 times.
✓ Branch 30 taken 55000 times.
✓ Branch 31 taken 118000 times.
✓ Branch 32 taken 55000 times.
✓ Branch 33 taken 118000 times.
✓ Branch 34 taken 55000 times.
✓ Branch 35 taken 118000 times.
✓ Branch 36 taken 55000 times.
✓ Branch 37 taken 118000 times.
✓ Branch 38 taken 55000 times.
✓ Branch 39 taken 118000 times.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
|
2362160 | if (globalPos[2] > this->gridGeometry().bBoxMax()[2] - eps_ || globalPos[2] < this->gridGeometry().bBoxMin()[2] + eps_) |
| 104 | values.setAllNeumann(); | ||
| 105 | else | ||
| 106 | values.setAllDirichlet(); | ||
| 107 | return values; | ||
| 108 | } | ||
| 109 | |||
| 110 | /*! | ||
| 111 | * \brief Evaluates the boundary conditions for a Dirichlet boundary segment. | ||
| 112 | * | ||
| 113 | * \param globalPos The position for which the bc type should be evaluated | ||
| 114 | * | ||
| 115 | * For this method, the \a values parameter stores primary variables. | ||
| 116 | */ | ||
| 117 | PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const | ||
| 118 | { | ||
| 119 | 128656 | PrimaryVariables values; | |
| 120 | 128656 | values = exactSolution(globalPos); | |
| 121 | return values; | ||
| 122 | } | ||
| 123 | |||
| 124 | /*! | ||
| 125 | * \brief Evaluate the boundary conditions for a neumann | ||
| 126 | * boundary segment. | ||
| 127 | * | ||
| 128 | * This is the method for the case where the Neumann condition is | ||
| 129 | * potentially solution dependent | ||
| 130 | * | ||
| 131 | * \param element The finite element | ||
| 132 | * \param fvGeometry The finite-volume geometry | ||
| 133 | * \param elemVolVars All volume variables for the element | ||
| 134 | * \param scvf The sub control volume face | ||
| 135 | * | ||
| 136 | * Negative values mean influx. | ||
| 137 | * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$. | ||
| 138 | */ | ||
| 139 | template<class ElementVolumeVariables, class ElemFluxVarsCache> | ||
| 140 | 167456 | NumEqVector neumann(const Element& element, | |
| 141 | const FVElementGeometry& fvGeometry, | ||
| 142 | const ElementVolumeVariables& elemVolVars, | ||
| 143 | const ElemFluxVarsCache&, | ||
| 144 | const SubControlVolumeFace& scvf) const | ||
| 145 | { | ||
| 146 | 167456 | NumEqVector flux(0.0); | |
| 147 | // integrate over the scvf to compute the flux | ||
| 148 | 167456 | const auto geometry = fvGeometry.geometry(scvf); | |
| 149 | 167456 | Scalar derivative = 0.0; | |
| 150 | 334912 | const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension-1>::rule(geometry.type(), 4); | |
| 151 |
4/4✓ Branch 0 taken 1507104 times.
✓ Branch 1 taken 167456 times.
✓ Branch 2 taken 1507104 times.
✓ Branch 3 taken 167456 times.
|
3516576 | for(auto&& qp : quad) |
| 152 | { | ||
| 153 | 3014208 | auto globalPos = geometry.global(qp.position()); | |
| 154 | 1507104 | globalPos[2] = 0; // the derivative in z-direction is the exact solution evaluated with z=0 | |
| 155 | 2497104 | derivative += exactSolution(globalPos)*qp.weight()*geometry.integrationElement(qp.position()); | |
| 156 | } | ||
| 157 | |||
| 158 | 167456 | const auto globalPos = scvf.ipGlobal(); | |
| 159 |
10/10✓ Branch 0 taken 80992 times.
✓ Branch 1 taken 86464 times.
✓ Branch 2 taken 80992 times.
✓ Branch 3 taken 86464 times.
✓ Branch 4 taken 80992 times.
✓ Branch 5 taken 86464 times.
✓ Branch 6 taken 80992 times.
✓ Branch 7 taken 86464 times.
✓ Branch 8 taken 80992 times.
✓ Branch 9 taken 86464 times.
|
837280 | if (globalPos[2] > this->gridGeometry().bBoxMax()[2] - eps_ ) |
| 160 | 80992 | flux[0] = -derivative/scvf.area(); | |
| 161 | else | ||
| 162 | 86464 | flux[0] = derivative/scvf.area(); | |
| 163 | 167456 | return flux; | |
| 164 | } | ||
| 165 | |||
| 166 | |||
| 167 | // \} | ||
| 168 | |||
| 169 | /*! | ||
| 170 | * \name Volume terms | ||
| 171 | */ | ||
| 172 | // \{ | ||
| 173 | |||
| 174 | /*! | ||
| 175 | * \brief Applies a vector of point sources which are possibly solution dependent. | ||
| 176 | * | ||
| 177 | * \param pointSources A vector of Dumux::PointSource s that contain | ||
| 178 | source values for all phases and space positions. | ||
| 179 | * | ||
| 180 | * For this method, the \a values method of the point source | ||
| 181 | * has to return the absolute mass rate in kg/s. Positive values mean | ||
| 182 | * that mass is created, negative ones mean that it vanishes. | ||
| 183 | */ | ||
| 184 | void addPointSources(std::vector<PointSource>& pointSources) const | ||
| 185 |
1/2✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
|
19 | { pointSources = this->couplingManager().bulkPointSources(); } |
| 186 | |||
| 187 | /*! | ||
| 188 | * \brief Evaluates the point sources (added by addPointSources) | ||
| 189 | * for all phases within a given sub control volume. | ||
| 190 | * | ||
| 191 | * This is the method for the case where the point source is | ||
| 192 | * solution dependent and requires some quantities that | ||
| 193 | * are specific to the fully-implicit method. | ||
| 194 | * | ||
| 195 | * \param source A single point source | ||
| 196 | * \param element The finite element | ||
| 197 | * \param fvGeometry The finite-volume geometry | ||
| 198 | * \param elemVolVars All volume variables for the element | ||
| 199 | * \param scv The sub-control volume within the element | ||
| 200 | * | ||
| 201 | * For this method, the \a values() method of the point sources returns | ||
| 202 | * the absolute rate mass generated or annihilated in kg/s. Positive values mean | ||
| 203 | * that mass is created, negative ones mean that it vanishes. | ||
| 204 | */ | ||
| 205 | template<class ElementVolumeVariables> | ||
| 206 | 1267416 | void pointSource(PointSource& source, | |
| 207 | const Element &element, | ||
| 208 | const FVElementGeometry& fvGeometry, | ||
| 209 | const ElementVolumeVariables& elemVolVars, | ||
| 210 | const SubControlVolume &scv) const | ||
| 211 | { | ||
| 212 | if constexpr (CouplingManager::couplingMode != Embedded1d3dCouplingMode::kernel) | ||
| 213 | { | ||
| 214 | // compute source at every integration point | ||
| 215 | 3802248 | const Scalar pressure3D = this->couplingManager().bulkPriVars(source.id())[Indices::pressureIdx]; | |
| 216 | 3802248 | const Scalar pressure1D = this->couplingManager().lowDimPriVars(source.id())[Indices::pressureIdx]; | |
| 217 | |||
| 218 | // calculate the source | ||
| 219 | 3802248 | const Scalar radius = this->couplingManager().radius(source.id()); | |
| 220 | 1267416 | const Scalar beta = 2*M_PI/(2*M_PI + std::log(radius)); | |
| 221 | 1267416 | const Scalar sourceValue = beta*(pressure1D - pressure3D); | |
| 222 | 2534832 | source = sourceValue*source.quadratureWeight()*source.integrationElement(); | |
| 223 | } | ||
| 224 | 1267416 | } | |
| 225 | |||
| 226 | /*! | ||
| 227 | * \brief Evaluates the source term for all phases within a given | ||
| 228 | * sub control volume. | ||
| 229 | * | ||
| 230 | * This is the method for the case where the source term is | ||
| 231 | * potentially solution dependent and requires some quantities that | ||
| 232 | * are specific to the fully-implicit method. | ||
| 233 | * | ||
| 234 | * \param element The finite element | ||
| 235 | * \param fvGeometry The finite-volume geometry | ||
| 236 | * \param elemVolVars All volume variables for the element | ||
| 237 | * \param scv The sub control volume | ||
| 238 | * | ||
| 239 | * For this method, the return parameter stores the conserved quantity rate | ||
| 240 | * generated or annihilated per volume unit. Positive values mean | ||
| 241 | * that the conserved quantity is created, negative ones mean that it vanishes. | ||
| 242 | * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$. | ||
| 243 | */ | ||
| 244 | template<class ElementVolumeVariables> | ||
| 245 | 925620 | NumEqVector source(const Element &element, | |
| 246 | const FVElementGeometry& fvGeometry, | ||
| 247 | const ElementVolumeVariables& elemVolVars, | ||
| 248 | const SubControlVolume &scv) const | ||
| 249 | { | ||
| 250 |
0/7✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
5544308 | NumEqVector source(0.0); |
| 251 | |||
| 252 | if constexpr(CouplingManager::couplingMode == Embedded1d3dCouplingMode::kernel) | ||
| 253 | { | ||
| 254 | 2776860 | const auto eIdx = this->gridGeometry().elementMapper().index(element); | |
| 255 | 2776860 | const auto& sourceIds = this->couplingManager().bulkSourceIds(eIdx, scv.indexInElement()); | |
| 256 | 2776860 | const auto& sourceWeights = this->couplingManager().bulkSourceWeights(eIdx, scv.indexInElement()); | |
| 257 | |||
| 258 |
4/4✓ Branch 0 taken 138480 times.
✓ Branch 1 taken 925620 times.
✓ Branch 2 taken 138480 times.
✓ Branch 3 taken 925620 times.
|
1202580 | for (int i = 0; i < sourceIds.size(); ++i) |
| 259 | { | ||
| 260 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 138480 times.
|
138480 | const auto id = sourceIds[i]; |
| 261 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 138480 times.
|
138480 | const auto weight = sourceWeights[i]; |
| 262 | 415440 | const auto xi = this->couplingManager().fluxScalingFactor(id); | |
| 263 | |||
| 264 | 415440 | const Scalar radius = this->couplingManager().radius(id); | |
| 265 | 415440 | const Scalar pressure3D = this->couplingManager().bulkPriVars(id)[Indices::pressureIdx]; | |
| 266 | 415440 | const Scalar pressure1D = this->couplingManager().lowDimPriVars(id)[Indices::pressureIdx]; | |
| 267 | |||
| 268 | // calculate the source | ||
| 269 |
3/4✓ Branch 0 taken 5 times.
✓ Branch 1 taken 138475 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
|
138480 | static const Scalar beta = 2*M_PI/(2*M_PI + std::log(radius)); |
| 270 | 138480 | const Scalar sourceValue = beta*(pressure1D - pressure3D)*xi; | |
| 271 | 138480 | source[Indices::conti0EqIdx] += sourceValue*weight; | |
| 272 | } | ||
| 273 | |||
| 274 | 925620 | const auto volume = scv.volume()*elemVolVars[scv].extrusionFactor(); | |
| 275 | 925620 | source[Indices::conti0EqIdx] /= volume; | |
| 276 | } | ||
| 277 | |||
| 278 | 925620 | return source; | |
| 279 | } | ||
| 280 | |||
| 281 | //! compute the flux scaling factor (xi) for a distance r for a vessel with radius R and kernel width rho | ||
| 282 | 380 | Scalar fluxScalingFactor(const Scalar r, const Scalar R, const Scalar rho) const | |
| 283 | { | ||
| 284 | using std::log; | ||
| 285 |
3/4✓ Branch 0 taken 5 times.
✓ Branch 1 taken 375 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
|
380 | static const Scalar beta = 2*M_PI/(2*M_PI + log(R)); |
| 286 |
3/4✓ Branch 0 taken 5 times.
✓ Branch 1 taken 375 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
|
380 | static const Scalar kernelWidthFactorLog = log(rho/R); |
| 287 |
2/2✓ Branch 0 taken 360 times.
✓ Branch 1 taken 20 times.
|
380 | return r < rho ? 1.0/(1.0 + R*beta/(2*M_PI*R)*(r*r/(2*rho*rho) + kernelWidthFactorLog - 0.5)) |
| 288 | 20 | : 1.0/(1.0 + R*beta/(2*M_PI*R)*log(r/R)); | |
| 289 | } | ||
| 290 | |||
| 291 | /*! | ||
| 292 | * \brief Evaluates the initial value for a control volume. | ||
| 293 | * | ||
| 294 | * \param globalPos The position for which the initial condition should be evaluated | ||
| 295 | * | ||
| 296 | * For this method, the \a values parameter stores primary | ||
| 297 | * variables. | ||
| 298 | */ | ||
| 299 | ✗ | PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const | |
| 300 | 933500 | { return PrimaryVariables(0.0); } | |
| 301 | |||
| 302 | ✗ | Scalar p1DExact(const GlobalPosition &globalPos) const | |
| 303 | 7188444 | { return 1.0 + globalPos[2]; } | |
| 304 | |||
| 305 | //! The exact solution | ||
| 306 | 3594222 | Scalar exactSolution(const GlobalPosition &globalPos) const | |
| 307 | { | ||
| 308 |
4/4✓ Branch 0 taken 19 times.
✓ Branch 1 taken 3594203 times.
✓ Branch 2 taken 19 times.
✓ Branch 3 taken 3594203 times.
|
7188444 | const auto r = std::hypot(globalPos[0], globalPos[1]); |
| 309 |
4/6✓ Branch 0 taken 19 times.
✓ Branch 1 taken 3594203 times.
✓ Branch 3 taken 19 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 19 times.
✗ Branch 7 not taken.
|
3594222 | static const auto R = getParam<Scalar>("SpatialParams.Radius"); |
| 310 | |||
| 311 | if (CouplingManager::couplingMode == Embedded1d3dCouplingMode::kernel) | ||
| 312 | { | ||
| 313 |
4/6✓ Branch 0 taken 5 times.
✓ Branch 1 taken 958595 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
|
958600 | static const auto rho = getParam<Scalar>("MixedDimension.KernelWidthFactor")*R; |
| 314 |
2/2✓ Branch 0 taken 949208 times.
✓ Branch 1 taken 9392 times.
|
958600 | if (r > rho) |
| 315 | 1898416 | return -1.0*p1DExact(globalPos)/(2*M_PI)*std::log(r); | |
| 316 | else | ||
| 317 | 18784 | return -1.0*p1DExact(globalPos)/(2*M_PI)*(r*r/(2.0*rho*rho) + std::log(rho) - 0.5); | |
| 318 | } | ||
| 319 | else if (CouplingManager::couplingMode == Embedded1d3dCouplingMode::surface) | ||
| 320 | { | ||
| 321 |
2/2✓ Branch 0 taken 949208 times.
✓ Branch 1 taken 9392 times.
|
958600 | if (r > R) |
| 322 | 1898416 | return -1.0*p1DExact(globalPos)/(2*M_PI)*std::log(r); | |
| 323 | else | ||
| 324 | 18784 | return -1.0*p1DExact(globalPos)/(2*M_PI)*std::log(R); | |
| 325 | } | ||
| 326 | |||
| 327 | 3354044 | return -1.0*p1DExact(globalPos)/(2*M_PI)*std::log(r); | |
| 328 | } | ||
| 329 | |||
| 330 | //! Called after every time step | ||
| 331 | //! Output the total global exchange term | ||
| 332 | 19 | void computeSourceIntegral(const SolutionVector& sol, const GridVariables& gridVars) | |
| 333 | { | ||
| 334 | 19 | PrimaryVariables source(0.0); | |
| 335 |
1/2✓ Branch 3 taken 931237 times.
✗ Branch 4 not taken.
|
1862493 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 336 | { | ||
| 337 | 2793654 | const auto fvGeometry = localView(this->gridGeometry()).bindElement(element); | |
| 338 | 4628654 | const auto elemVolVars = localView(gridVars.curGridVolVars()).bindElement(element, fvGeometry, sol); | |
| 339 | |||
| 340 |
5/6✓ Branch 0 taken 1027244 times.
✓ Branch 1 taken 931218 times.
✓ Branch 2 taken 1027244 times.
✓ Branch 3 taken 931218 times.
✓ Branch 5 taken 917500 times.
✗ Branch 6 not taken.
|
3916924 | for (auto&& scv : scvs(fvGeometry)) |
| 341 | { | ||
| 342 |
1/2✓ Branch 1 taken 917500 times.
✗ Branch 2 not taken.
|
1027244 | auto pointSources = this->scvPointSources(element, fvGeometry, elemVolVars, scv); |
| 343 | 1753988 | pointSources += this->source(element, fvGeometry, elemVolVars, scv); | |
| 344 | 1136988 | pointSources *= scv.volume()*elemVolVars[scv].extrusionFactor(); | |
| 345 | 1027244 | source += pointSources; | |
| 346 | } | ||
| 347 | } | ||
| 348 | |||
| 349 | 38 | std::cout << "Global integrated source (3D): " << source << '\n'; | |
| 350 | 19 | } | |
| 351 | |||
| 352 | /*! | ||
| 353 | * \brief Adds additional VTK output data to the VTKWriter. | ||
| 354 | * | ||
| 355 | * Function is called by the output module on every write. | ||
| 356 | */ | ||
| 357 | template<class VtkOutputModule> | ||
| 358 | 19 | void addVtkOutputFields(VtkOutputModule& vtk) const | |
| 359 | { | ||
| 360 |
4/10✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 19 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 19 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
38 | vtk.addField(exactPressure_, "exact pressure"); |
| 361 | 19 | } | |
| 362 | |||
| 363 | //! Get the coupling manager | ||
| 364 | const CouplingManager& couplingManager() const | ||
| 365 |
15/21✗ Branch 0 not taken.
✓ Branch 1 taken 139280 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1367064 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1405896 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1405896 times.
✓ Branch 8 taken 5 times.
✓ Branch 9 taken 1367061 times.
✓ Branch 10 taken 5 times.
✓ Branch 11 taken 139286 times.
✓ Branch 12 taken 7 times.
✓ Branch 13 taken 138476 times.
✓ Branch 14 taken 16 times.
✓ Branch 15 taken 138475 times.
✓ Branch 16 taken 1 times.
✓ Branch 17 taken 5 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 5 times.
✗ Branch 21 not taken.
|
6794470 | { return *couplingManager_; } |
| 366 | |||
| 367 | private: | ||
| 368 | std::vector<Scalar> exactPressure_; | ||
| 369 | |||
| 370 | static constexpr Scalar eps_ = 1.5e-7; | ||
| 371 | std::string name_; | ||
| 372 | |||
| 373 | std::shared_ptr<CouplingManager> couplingManager_; | ||
| 374 | }; | ||
| 375 | |||
| 376 | } // end namespace Dumux | ||
| 377 | |||
| 378 | #endif | ||
| 379 |