| 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 Experimental | ||
| 10 | * \ingroup Assembly | ||
| 11 | * \brief A linear system assembler (residual and Jacobian) for finite volume schemes | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_EXPERIMENTAL_MULTISTAGE_FV_ASSEMBLER_HH | ||
| 14 | #define DUMUX_EXPERIMENTAL_MULTISTAGE_FV_ASSEMBLER_HH | ||
| 15 | |||
| 16 | #include <vector> | ||
| 17 | #include <deque> | ||
| 18 | #include <type_traits> | ||
| 19 | #include <memory> | ||
| 20 | |||
| 21 | #include <dune/istl/matrixindexset.hh> | ||
| 22 | |||
| 23 | #include <dumux/common/properties.hh> | ||
| 24 | #include <dumux/common/gridcapabilities.hh> | ||
| 25 | #include <dumux/common/typetraits/vector.hh> | ||
| 26 | |||
| 27 | #include <dumux/discretization/method.hh> | ||
| 28 | #include <dumux/linear/parallelhelpers.hh> | ||
| 29 | #include <dumux/linear/dunevectors.hh> | ||
| 30 | |||
| 31 | #include <dumux/assembly/coloring.hh> | ||
| 32 | #include <dumux/assembly/jacobianpattern.hh> | ||
| 33 | #include <dumux/assembly/diffmethod.hh> | ||
| 34 | |||
| 35 | #include <dumux/parallel/multithreading.hh> | ||
| 36 | #include <dumux/parallel/parallel_for.hh> | ||
| 37 | |||
| 38 | #include <dumux/experimental/assembly/cvfelocalassembler.hh> | ||
| 39 | #include <dumux/experimental/assembly/cclocalassembler.hh> | ||
| 40 | #include <dumux/experimental/assembly/multistagefvlocaloperator.hh> | ||
| 41 | |||
| 42 | #include <dumux/experimental/timestepping/multistagemethods.hh> | ||
| 43 | #include <dumux/experimental/timestepping/multistagetimestepper.hh> | ||
| 44 | |||
| 45 | namespace Dumux::Experimental::Detail { | ||
| 46 | |||
| 47 | template<class DiscretizationMethod> | ||
| 48 | struct LocalAssemblerChooser; | ||
| 49 | |||
| 50 | template<class DM> | ||
| 51 | struct LocalAssemblerChooser<DiscretizationMethods::CVFE<DM>> | ||
| 52 | { | ||
| 53 | template<class TypeTag, class Impl, DiffMethod diffMethod> | ||
| 54 | using type = Experimental::CVFELocalAssembler<TypeTag, Impl, diffMethod>; | ||
| 55 | }; | ||
| 56 | |||
| 57 | template<> | ||
| 58 | struct LocalAssemblerChooser<DiscretizationMethods::CCMpfa> | ||
| 59 | { | ||
| 60 | template<class TypeTag, class Impl, DiffMethod diffMethod> | ||
| 61 | using type = Experimental::CCLocalAssembler<TypeTag, Impl, diffMethod>; | ||
| 62 | }; | ||
| 63 | |||
| 64 | template<> | ||
| 65 | struct LocalAssemblerChooser<DiscretizationMethods::CCTpfa> | ||
| 66 | { | ||
| 67 | template<class TypeTag, class Impl, DiffMethod diffMethod> | ||
| 68 | using type = Experimental::CCLocalAssembler<TypeTag, Impl, diffMethod>; | ||
| 69 | }; | ||
| 70 | |||
| 71 | template<class TypeTag, class Impl, DiffMethod diffMethod> | ||
| 72 | using LocalAssemblerChooser_t = typename LocalAssemblerChooser< | ||
| 73 | typename GetPropType<TypeTag, Properties::GridGeometry>::DiscretizationMethod | ||
| 74 | >::template type<TypeTag, Impl, diffMethod>; | ||
| 75 | |||
| 76 | } // end namespace Dumux::Detail | ||
| 77 | |||
| 78 | namespace Dumux::Experimental { | ||
| 79 | |||
| 80 | /*! | ||
| 81 | * \ingroup Experimental | ||
| 82 | * \ingroup Assembly | ||
| 83 | * \brief A linear system assembler (residual and Jacobian) for finite volume schemes (box, tpfa, mpfa, ...) | ||
| 84 | * \tparam TypeTag The TypeTag | ||
| 85 | * \tparam diffMethod The differentiation method to residual compute derivatives | ||
| 86 | */ | ||
| 87 | template<class TypeTag, DiffMethod diffMethod> | ||
| 88 | class MultiStageFVAssembler | ||
| 89 | { | ||
| 90 | using GridGeo = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 91 | using GridView = typename GridGeo::GridView; | ||
| 92 | using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>; | ||
| 93 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 94 | using ElementSeed = typename GridView::Grid::template Codim<0>::EntitySeed; | ||
| 95 | |||
| 96 | static constexpr bool isBox = GridGeo::discMethod == DiscretizationMethods::box; | ||
| 97 | |||
| 98 | using ThisType = MultiStageFVAssembler<TypeTag, diffMethod>; | ||
| 99 | using LocalAssembler = typename Detail::LocalAssemblerChooser_t<TypeTag, ThisType, diffMethod>; | ||
| 100 | |||
| 101 | public: | ||
| 102 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 103 | using StageParams = Experimental::MultiStageParams<Scalar>; | ||
| 104 | using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>; | ||
| 105 | using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; | ||
| 106 | using ResidualType = typename Dumux::Detail::NativeDuneVectorType<SolutionVector>::type; | ||
| 107 | |||
| 108 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 109 | |||
| 110 | using GridGeometry = GridGeo; | ||
| 111 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 112 | |||
| 113 | /*! | ||
| 114 | * \brief The constructor for instationary problems | ||
| 115 | * \note the grid variables might be temporarily changed during assembly (if caching is enabled) | ||
| 116 | * it is however guaranteed that the state after assembly will be the same as before | ||
| 117 | */ | ||
| 118 | 9 | MultiStageFVAssembler(std::shared_ptr<const Problem> problem, | |
| 119 | std::shared_ptr<const GridGeometry> gridGeometry, | ||
| 120 | std::shared_ptr<GridVariables> gridVariables, | ||
| 121 | std::shared_ptr<const Experimental::MultiStageMethod<Scalar>> msMethod, | ||
| 122 | const SolutionVector& prevSol) | ||
| 123 | : timeSteppingMethod_(msMethod) | ||
| 124 | , problem_(problem) | ||
| 125 | , gridGeometry_(gridGeometry) | ||
| 126 | , gridVariables_(gridVariables) | ||
| 127 |
2/20✓ Branch 9 taken 9 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 9 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
|
9 | , prevSol_(&prevSol) |
| 128 | { | ||
| 129 | 9 | enableMultithreading_ = SupportsColoring<typename GridGeometry::DiscretizationMethod>::value | |
| 130 |
3/6✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
|
27 | && Grid::Capabilities::supportsMultithreading(gridGeometry_->gridView()) |
| 131 | && !Multithreading::isSerial() | ||
| 132 |
2/4✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9 times.
|
9 | && getParam<bool>("Assembly.Multithreading", true); |
| 133 | |||
| 134 |
1/2✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
|
9 | maybeComputeColors_(); |
| 135 | 9 | } | |
| 136 | |||
| 137 | /*! | ||
| 138 | * \brief Assembles the global Jacobian of the residual | ||
| 139 | * and the residual for the current solution. | ||
| 140 | */ | ||
| 141 | template<class PartialReassembler = DefaultPartialReassembler> | ||
| 142 | 2172 | void assembleJacobianAndResidual(const SolutionVector& curSol, const PartialReassembler* partialReassembler = nullptr) | |
| 143 | { | ||
| 144 | 2172 | resetJacobian_(partialReassembler); | |
| 145 | |||
| 146 | 2172 | resetResidual_(); | |
| 147 | 4344 | spatialOperatorEvaluations_.back() = 0.0; | |
| 148 | 4344 | temporalOperatorEvaluations_.back() = 0.0; | |
| 149 | |||
| 150 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 2172 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2172 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2172 times.
|
6516 | if (stageParams_->size() != spatialOperatorEvaluations_.size()) |
| 151 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Wrong number of residuals"); | |
| 152 | |||
| 153 | 2172 | assemble_([&](const Element& element) | |
| 154 | { | ||
| 155 |
1/4✓ Branch 1 taken 10368 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
12652448 | LocalAssembler localAssembler(*this, element, curSol); |
| 156 |
1/2✓ Branch 1 taken 5066048 times.
✗ Branch 2 not taken.
|
5066048 | localAssembler.assembleJacobianAndResidual( |
| 157 |
3/6✓ Branch 1 taken 5066048 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5066048 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5066048 times.
✗ Branch 8 not taken.
|
15198144 | *jacobian_, *residual_, *gridVariables_, |
| 158 |
2/4✓ Branch 1 taken 5066048 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5066048 times.
✗ Branch 5 not taken.
|
10132096 | *stageParams_, |
| 159 | temporalOperatorEvaluations_.back(), | ||
| 160 | spatialOperatorEvaluations_.back(), | ||
| 161 | 5066048 | constrainedDofs_, | |
| 162 |
1/2✓ Branch 1 taken 5066048 times.
✗ Branch 2 not taken.
|
5066048 | partialReassembler |
| 163 | ); | ||
| 164 | }); | ||
| 165 | |||
| 166 | // assemble the full residual for the time integration stage | ||
| 167 |
1/2✓ Branch 2 taken 2172 times.
✗ Branch 3 not taken.
|
6516 | auto constantResidualComponent = (*residual_); |
| 168 | constantResidualComponent = 0.0; | ||
| 169 |
4/4✓ Branch 0 taken 2172 times.
✓ Branch 1 taken 2172 times.
✓ Branch 2 taken 2172 times.
✓ Branch 3 taken 2172 times.
|
6516 | for (std::size_t k = 0; k < stageParams_->size()-1; ++k) |
| 170 | { | ||
| 171 |
3/6✓ Branch 0 taken 2172 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2172 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2172 times.
✗ Branch 5 not taken.
|
6516 | if (!stageParams_->skipTemporal(k)) |
| 172 | 8688 | constantResidualComponent.axpy(stageParams_->temporalWeight(k), temporalOperatorEvaluations_[k]); | |
| 173 |
6/6✓ Branch 0 taken 2000 times.
✓ Branch 1 taken 172 times.
✓ Branch 2 taken 2000 times.
✓ Branch 3 taken 172 times.
✓ Branch 4 taken 2000 times.
✓ Branch 5 taken 172 times.
|
6516 | if (!stageParams_->skipSpatial(k)) |
| 174 | 8000 | constantResidualComponent.axpy(stageParams_->spatialWeight(k), spatialOperatorEvaluations_[k]); | |
| 175 | } | ||
| 176 | |||
| 177 | // masked summation of constant residual component onto this stage's resiudal component | ||
| 178 |
2/2✓ Branch 0 taken 5068836 times.
✓ Branch 1 taken 2172 times.
|
5071008 | for (std::size_t i = 0; i < constantResidualComponent.size(); ++i) |
| 179 |
4/4✓ Branch 0 taken 10137672 times.
✓ Branch 1 taken 5068836 times.
✓ Branch 2 taken 10137672 times.
✓ Branch 3 taken 5068836 times.
|
30413016 | for (std::size_t ii = 0; ii < constantResidualComponent[i].size(); ++ii) |
| 180 |
6/6✓ Branch 0 taken 10135360 times.
✓ Branch 1 taken 2312 times.
✓ Branch 2 taken 10135360 times.
✓ Branch 3 taken 2312 times.
✓ Branch 4 taken 10135360 times.
✓ Branch 5 taken 2312 times.
|
30413016 | (*residual_)[i][ii] += constrainedDofs_[i][ii] > 0.5 ? 0.0 : constantResidualComponent[i][ii]; |
| 181 | 2172 | } | |
| 182 | |||
| 183 | //! compute the residuals using the internal residual | ||
| 184 | ✗ | void assembleResidual(const SolutionVector& curSol) | |
| 185 | ✗ | { DUNE_THROW(Dune::NotImplemented, "residual"); } | |
| 186 | |||
| 187 | /*! | ||
| 188 | * \brief The version without arguments uses the default constructor to create | ||
| 189 | * the jacobian and residual objects in this assembler if you don't need them outside this class | ||
| 190 | */ | ||
| 191 | 9 | void setLinearSystem() | |
| 192 | { | ||
| 193 |
2/4✗ Branch 1 not taken.
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 9 times.
|
9 | jacobian_ = std::make_shared<JacobianMatrix>(); |
| 194 | 18 | jacobian_->setBuildMode(JacobianMatrix::random); | |
| 195 |
2/4✗ Branch 1 not taken.
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 9 times.
|
9 | residual_ = std::make_shared<ResidualType>(); |
| 196 | |||
| 197 | 15 | setResidualSize_(*residual_); | |
| 198 | 9 | setJacobianPattern_(); | |
| 199 | 9 | } | |
| 200 | |||
| 201 | /*! | ||
| 202 | * \brief Resizes jacobian and residual and recomputes colors | ||
| 203 | */ | ||
| 204 | void updateAfterGridAdaption() | ||
| 205 | { | ||
| 206 | setResidualSize_(*residual_); | ||
| 207 | setJacobianPattern_(); | ||
| 208 | maybeComputeColors_(); | ||
| 209 | } | ||
| 210 | |||
| 211 | //! Returns the number of degrees of freedom | ||
| 212 | std::size_t numDofs() const | ||
| 213 | 6294 | { return gridGeometry_->numDofs(); } | |
| 214 | |||
| 215 | //! The problem | ||
| 216 | const Problem& problem() const | ||
| 217 |
14/33✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 4 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 4 times.
✓ Branch 21 taken 3 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 4 times.
✓ Branch 24 taken 3 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 4 times.
✓ Branch 27 taken 3 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 3 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 3 times.
✗ Branch 34 not taken.
✓ Branch 36 taken 3 times.
✗ Branch 37 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
|
23790698 | { return *problem_; } |
| 218 | |||
| 219 | //! The global finite volume geometry | ||
| 220 | const GridGeometry& gridGeometry() const | ||
| 221 |
6/17✓ Branch 1 taken 10000000 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 10000000 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 52224 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 52224 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 20352 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 20352 times.
✗ Branch 19 not taken.
|
20145152 | { return *gridGeometry_; } |
| 222 | |||
| 223 | //! The gridview | ||
| 224 | const GridView& gridView() const | ||
| 225 |
4/16✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2172 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 5066048 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 2052 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 5019968 times.
|
10090240 | { return gridGeometry().gridView(); } |
| 226 | |||
| 227 | //! The global grid variables | ||
| 228 | GridVariables& gridVariables() | ||
| 229 | 4448 | { return *gridVariables_; } | |
| 230 | |||
| 231 | //! The global grid variables | ||
| 232 | const GridVariables& gridVariables() const | ||
| 233 |
4/8✓ Branch 1 taken 10072576 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10072576 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10072576 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 10072576 times.
✗ Branch 11 not taken.
|
40344064 | { return *gridVariables_; } |
| 234 | |||
| 235 | //! The jacobian matrix | ||
| 236 | JacobianMatrix& jacobian() | ||
| 237 | 8688 | { return *jacobian_; } | |
| 238 | |||
| 239 | //! The residual vector (rhs) | ||
| 240 | ResidualType& residual() | ||
| 241 |
2/4✓ Branch 3 taken 52 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 52 times.
✗ Branch 7 not taken.
|
4448 | { return *residual_; } |
| 242 | |||
| 243 | //! The solution of the previous time step | ||
| 244 | ✗ | const SolutionVector& prevSol() const | |
| 245 | ✗ | { return *prevSol_; } | |
| 246 | |||
| 247 | /*! | ||
| 248 | * \brief Create a local residual object (used by the local assembler) | ||
| 249 | */ | ||
| 250 | MultiStageFVLocalOperator<LocalResidual> localResidual() const | ||
| 251 |
3/6✓ Branch 1 taken 10072576 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10072576 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5059520 times.
✗ Branch 8 not taken.
|
25231552 | { return { LocalResidual{problem_.get(), nullptr} }; } |
| 252 | |||
| 253 | /*! | ||
| 254 | * \brief Update the grid variables | ||
| 255 | */ | ||
| 256 | ✗ | void updateGridVariables(const SolutionVector &cursol) | |
| 257 |
2/4✓ Branch 1 taken 2000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2000 times.
✗ Branch 5 not taken.
|
4158 | { gridVariables().update(cursol); } |
| 258 | |||
| 259 | /*! | ||
| 260 | * \brief Reset the gridVariables | ||
| 261 | */ | ||
| 262 | ✗ | void resetTimeStep(const SolutionVector &cursol) | |
| 263 | ✗ | { gridVariables().resetTimeStep(cursol); } | |
| 264 | |||
| 265 | 4104 | void clearStages() | |
| 266 | { | ||
| 267 | 4104 | spatialOperatorEvaluations_.clear(); | |
| 268 | 4104 | temporalOperatorEvaluations_.clear(); | |
| 269 |
2/2✓ Branch 0 taken 2052 times.
✓ Branch 1 taken 2052 times.
|
4104 | stageParams_.reset(); |
| 270 | 4104 | } | |
| 271 | |||
| 272 | template<class StageParams> | ||
| 273 | 2052 | void prepareStage(SolutionVector& x, StageParams params) | |
| 274 | { | ||
| 275 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2052 times.
|
2052 | stageParams_ = std::move(params); |
| 276 |
1/2✓ Branch 0 taken 2052 times.
✗ Branch 1 not taken.
|
2052 | const auto curStage = stageParams_->size() - 1; |
| 277 | |||
| 278 | // in the first stage, also assemble the residual | ||
| 279 | // at the previous time level (stage 0 residual) | ||
| 280 |
1/2✓ Branch 0 taken 2052 times.
✗ Branch 1 not taken.
|
2052 | if (curStage == 1) |
| 281 | { | ||
| 282 | // update time in variables? | ||
| 283 | 6156 | setProblemTime_(*problem_, stageParams_->timeAtStage(0)); | |
| 284 | |||
| 285 | 2052 | resetResidual_(); // residual resized and zero | |
| 286 | |||
| 287 | 2052 | assert(spatialOperatorEvaluations_.size() >= 0); | |
| 288 |
1/2✓ Branch 0 taken 2052 times.
✗ Branch 1 not taken.
|
2052 | if (spatialOperatorEvaluations_.size() == 0) |
| 289 | { | ||
| 290 | 4104 | spatialOperatorEvaluations_.push_back(*residual_); | |
| 291 | 4104 | temporalOperatorEvaluations_.push_back(*residual_); | |
| 292 | |||
| 293 | // assemble stage 0 residuals | ||
| 294 | 20112644 | assemble_([&](const auto& element) | |
| 295 | { | ||
| 296 |
1/5✓ Branch 1 taken 3072 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
12546080 | LocalAssembler localAssembler(*this, element, *prevSol_); |
| 297 |
1/4✓ Branch 1 taken 5019968 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
5019968 | localAssembler.localResidual().spatialWeight(1.0); |
| 298 |
1/4✓ Branch 1 taken 5019968 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
5019968 | localAssembler.localResidual().temporalWeight(1.0); |
| 299 |
3/12✓ Branch 1 taken 5019968 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5019968 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5019968 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
|
15059904 | localAssembler.assembleCurrentResidual(spatialOperatorEvaluations_.back(), temporalOperatorEvaluations_.back()); |
| 300 | }); | ||
| 301 | } | ||
| 302 | |||
| 303 | // we don't delete the first stage so it can be reused in a restarted | ||
| 304 | // time integration step. The evaluations are only deleted | ||
| 305 | // when explicitly requested by calling clearStages(). | ||
| 306 | // So if here the vector is non-empty, we don't need to evaluate again | ||
| 307 | // (this should only occur if we are restarting time integration, e.g. | ||
| 308 | // with a different time step size) | ||
| 309 | ✗ | else if (spatialOperatorEvaluations_.size() > 0) | |
| 310 | { | ||
| 311 | ✗ | updateGridVariables(x); | |
| 312 | ✗ | spatialOperatorEvaluations_.resize(1); | |
| 313 | ✗ | temporalOperatorEvaluations_.resize(1); | |
| 314 | } | ||
| 315 | } | ||
| 316 | |||
| 317 | // update time in variables? | ||
| 318 | 6156 | setProblemTime_(*problem_, stageParams_->timeAtStage(curStage)); | |
| 319 | |||
| 320 | 2052 | resetResidual_(); // residual resized and zero | |
| 321 | |||
| 322 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2052 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2052 times.
|
4104 | if (spatialOperatorEvaluations_.size() != curStage) |
| 323 | ✗ | DUNE_THROW(Dune::InvalidStateException, | |
| 324 | "Invalid state. Maybe you forgot to call clearStages()"); | ||
| 325 | |||
| 326 | // allocate memory for this stage | ||
| 327 | 4104 | spatialOperatorEvaluations_.push_back(*residual_); | |
| 328 | 4104 | temporalOperatorEvaluations_.push_back(*residual_); | |
| 329 | 2052 | } | |
| 330 | |||
| 331 | //! TODO get rid of this (called by Newton but shouldn't be necessary) | ||
| 332 | ✗ | bool isStationaryProblem() const | |
| 333 | ✗ | { return false; } | |
| 334 | |||
| 335 | bool isImplicit() const | ||
| 336 | 20172032 | { return timeSteppingMethod_->implicit(); } | |
| 337 | |||
| 338 | private: | ||
| 339 | /*! | ||
| 340 | * \brief Resizes the jacobian and sets the jacobian' sparsity pattern. | ||
| 341 | */ | ||
| 342 | 9 | void setJacobianPattern_() | |
| 343 | { | ||
| 344 | // resize the jacobian and the residual | ||
| 345 | 9 | const auto numDofs = this->numDofs(); | |
| 346 | 18 | jacobian_->setSize(numDofs, numDofs); | |
| 347 | |||
| 348 | // create occupation pattern of the jacobian | ||
| 349 |
2/2✓ Branch 2 taken 7 times.
✓ Branch 3 taken 2 times.
|
18 | if (timeSteppingMethod_->implicit()) |
| 350 |
2/4✓ Branch 3 taken 7 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 7 times.
✗ Branch 7 not taken.
|
14 | getJacobianPattern<true>(gridGeometry()).exportIdx(*jacobian_); |
| 351 | else | ||
| 352 |
2/4✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
|
4 | getJacobianPattern<false>(gridGeometry()).exportIdx(*jacobian_); |
| 353 | 9 | } | |
| 354 | |||
| 355 | //! Resizes the residual | ||
| 356 | 6174 | void setResidualSize_(ResidualType& res) | |
| 357 | 6288 | { res.resize(numDofs()); } | |
| 358 | |||
| 359 | //! Computes the colors | ||
| 360 | 9 | void maybeComputeColors_() | |
| 361 | { | ||
| 362 |
1/2✓ Branch 0 taken 9 times.
✗ Branch 1 not taken.
|
9 | if (enableMultithreading_) |
| 363 | 18 | elementSets_ = computeColoring(gridGeometry()).sets; | |
| 364 | 9 | } | |
| 365 | |||
| 366 | // reset the residual vector to 0.0 | ||
| 367 | 6276 | void resetResidual_() | |
| 368 | { | ||
| 369 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6276 times.
|
6276 | if(!residual_) |
| 370 | { | ||
| 371 | ✗ | residual_ = std::make_shared<ResidualType>(); | |
| 372 | ✗ | setResidualSize_(*residual_); | |
| 373 | } | ||
| 374 | |||
| 375 | 6276 | setResidualSize_(constrainedDofs_); | |
| 376 | |||
| 377 | 12552 | (*residual_) = 0.0; | |
| 378 | 6276 | constrainedDofs_ = 0.0; | |
| 379 | 6276 | } | |
| 380 | |||
| 381 | // reset the Jacobian matrix to 0.0 | ||
| 382 | template <class PartialReassembler = DefaultPartialReassembler> | ||
| 383 | 2172 | void resetJacobian_(const PartialReassembler *partialReassembler = nullptr) | |
| 384 | { | ||
| 385 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2172 times.
|
2172 | if(!jacobian_) |
| 386 | { | ||
| 387 | ✗ | jacobian_ = std::make_shared<JacobianMatrix>(); | |
| 388 | ✗ | jacobian_->setBuildMode(JacobianMatrix::random); | |
| 389 | ✗ | setJacobianPattern_(); | |
| 390 | } | ||
| 391 | |||
| 392 |
2/2✓ Branch 0 taken 2000 times.
✓ Branch 1 taken 172 times.
|
2172 | if (partialReassembler) |
| 393 | ✗ | partialReassembler->resetJacobian(*this); | |
| 394 | else | ||
| 395 | 4344 | *jacobian_ = 0.0; | |
| 396 | 2172 | } | |
| 397 | |||
| 398 | /*! | ||
| 399 | * \brief A method assembling something per element | ||
| 400 | * \note Handles exceptions for parallel runs | ||
| 401 | * \throws NumericalProblem on all processes if an exception is thrown during assembly | ||
| 402 | */ | ||
| 403 | template<typename AssembleElementFunc> | ||
| 404 | 8448 | void assemble_(AssembleElementFunc&& assembleElement) const | |
| 405 | { | ||
| 406 | // a state that will be checked on all processes | ||
| 407 | 8448 | bool succeeded = false; | |
| 408 | |||
| 409 | // try assembling using the local assembly function | ||
| 410 | try | ||
| 411 | { | ||
| 412 |
1/2✓ Branch 0 taken 4224 times.
✗ Branch 1 not taken.
|
8448 | if (enableMultithreading_) |
| 413 | { | ||
| 414 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 4224 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4224 times.
|
16896 | assert(elementSets_.size() > 0); |
| 415 | |||
| 416 | // make this element loop run in parallel | ||
| 417 | // for this we have to color the elements so that we don't get | ||
| 418 | // race conditions when writing into the global matrix | ||
| 419 | // each color can be assembled using multiple threads | ||
| 420 |
4/4✓ Branch 0 taken 37508 times.
✓ Branch 1 taken 4224 times.
✓ Branch 2 taken 37508 times.
✓ Branch 3 taken 4224 times.
|
166928 | for (const auto& elements : elementSets_) |
| 421 | { | ||
| 422 |
1/2✓ Branch 2 taken 37508 times.
✗ Branch 3 not taken.
|
20322064 | Dumux::parallelFor(elements.size(), [&](const std::size_t i) |
| 423 | { | ||
| 424 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 5066048 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5066048 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 5019968 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 5019968 times.
|
20172032 | const auto element = gridView().grid().entity(elements[i]); |
| 425 | 10086016 | assembleElement(element); | |
| 426 | }); | ||
| 427 | } | ||
| 428 | } | ||
| 429 | else | ||
| 430 | ✗ | for (const auto& element : elements(gridView())) | |
| 431 | ✗ | assembleElement(element); | |
| 432 | |||
| 433 | // if we get here, everything worked well on this process | ||
| 434 | 8448 | succeeded = true; | |
| 435 | } | ||
| 436 | // throw exception if a problem occurred | ||
| 437 | ✗ | catch (NumericalProblem &e) | |
| 438 | { | ||
| 439 | ✗ | std::cout << "rank " << gridView().comm().rank() | |
| 440 | ✗ | << " caught an exception while assembling:" << e.what() | |
| 441 | ✗ | << "\n"; | |
| 442 | ✗ | succeeded = false; | |
| 443 | } | ||
| 444 | |||
| 445 | // make sure everything worked well on all processes | ||
| 446 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 4224 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4224 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4224 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4224 times.
|
16896 | if (gridView().comm().size() > 1) |
| 447 | ✗ | succeeded = gridView().comm().min(succeeded); | |
| 448 | |||
| 449 | // if not succeeded rethrow the error on all processes | ||
| 450 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4224 times.
|
8448 | if (!succeeded) |
| 451 | ✗ | DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system"); | |
| 452 | 8448 | } | |
| 453 | |||
| 454 | // TODO make this nicer with a is_detected trait in a common location | ||
| 455 | template<class P> | ||
| 456 | ✗ | void setProblemTime_(const P& p, const Scalar t) | |
| 457 | 4104 | { setProblemTimeImpl_(p, t, 0); } | |
| 458 | |||
| 459 | template<class P> | ||
| 460 | auto setProblemTimeImpl_(const P& p, const Scalar t, int) -> decltype(p.setTime(0)) | ||
| 461 | { p.setTime(t); } | ||
| 462 | |||
| 463 | template<class P> | ||
| 464 | ✗ | void setProblemTimeImpl_(const P& p, const Scalar t, long) | |
| 465 | ✗ | {} | |
| 466 | |||
| 467 | std::shared_ptr<const Experimental::MultiStageMethod<Scalar>> timeSteppingMethod_; | ||
| 468 | std::vector<ResidualType> spatialOperatorEvaluations_; | ||
| 469 | std::vector<ResidualType> temporalOperatorEvaluations_; | ||
| 470 | ResidualType constrainedDofs_; | ||
| 471 | std::shared_ptr<const StageParams> stageParams_; | ||
| 472 | |||
| 473 | //! pointer to the problem to be solved | ||
| 474 | std::shared_ptr<const Problem> problem_; | ||
| 475 | |||
| 476 | //! the finite volume geometry of the grid | ||
| 477 | std::shared_ptr<const GridGeometry> gridGeometry_; | ||
| 478 | |||
| 479 | //! the variables container for the grid | ||
| 480 | std::shared_ptr<GridVariables> gridVariables_; | ||
| 481 | |||
| 482 | //! an observing pointer to the previous solution for instationary problems | ||
| 483 | const SolutionVector* prevSol_ = nullptr; | ||
| 484 | |||
| 485 | //! shared pointers to the jacobian matrix and residual | ||
| 486 | std::shared_ptr<JacobianMatrix> jacobian_; | ||
| 487 | std::shared_ptr<ResidualType> residual_; | ||
| 488 | |||
| 489 | //! element sets for parallel assembly | ||
| 490 | bool enableMultithreading_ = false; | ||
| 491 | std::deque<std::vector<ElementSeed>> elementSets_; | ||
| 492 | }; | ||
| 493 | |||
| 494 | } // namespace Dumux | ||
| 495 | |||
| 496 | #endif | ||
| 497 |