| 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 | // In the file `main.cc`, we use the previously defined model to | ||
| 9 | // setup the simulation. The setup consists of four steps: | ||
| 10 | // 1. Define the problem setting boundary conditions and the diffusion coefficient | ||
| 11 | // 2. Configure the property system reusing the model defined in Part 1 | ||
| 12 | // 3. Define a function for setting the random initial condition | ||
| 13 | // 4. The main program defining all steps of the program | ||
| 14 | // | ||
| 15 | // __Table of contents__ | ||
| 16 | // | ||
| 17 | // [TOC] | ||
| 18 | // | ||
| 19 | // We start in `main.cc` with the necessary header includes: | ||
| 20 | // [[details]] includes | ||
| 21 | #include <config.h> | ||
| 22 | |||
| 23 | // Include the header containing std::true_type and std::false_type | ||
| 24 | #include <type_traits> | ||
| 25 | |||
| 26 | // Include the headers useful for the random initial solution | ||
| 27 | #include <random> | ||
| 28 | #include <dumux/common/random.hh> | ||
| 29 | |||
| 30 | // As Dune grid interface implementation we will use | ||
| 31 | // the structured parallel grid manager YaspGrid | ||
| 32 | #include <dune/grid/yaspgrid.hh> | ||
| 33 | |||
| 34 | // Common includes for problem and main | ||
| 35 | // [[codeblock]] | ||
| 36 | #include <dumux/common/initialize.hh> | ||
| 37 | #include <dumux/common/properties.hh> | ||
| 38 | #include <dumux/common/parameters.hh> | ||
| 39 | #include <dumux/common/numeqvector.hh> | ||
| 40 | #include <dumux/common/boundarytypes.hh> | ||
| 41 | #include <dumux/common/fvproblem.hh> | ||
| 42 | |||
| 43 | // VTK output functionality | ||
| 44 | #include <dumux/io/vtkoutputmodule.hh> | ||
| 45 | #include <dumux/io/grid/gridmanager_yasp.hh> | ||
| 46 | |||
| 47 | // the discretization method | ||
| 48 | #include <dumux/discretization/box.hh> | ||
| 49 | |||
| 50 | // assembly and solver | ||
| 51 | #include <dumux/linear/linearsolvertraits.hh> | ||
| 52 | #include <dumux/linear/linearalgebratraits.hh> | ||
| 53 | #include <dumux/linear/istlsolvers.hh> | ||
| 54 | #include <dumux/linear/pdesolver.hh> | ||
| 55 | #include <dumux/assembly/fvassembler.hh> | ||
| 56 | // [[/codeblock]] | ||
| 57 | // | ||
| 58 | // Finally, we include the model defined in Part 1. | ||
| 59 | #include "model.hh" | ||
| 60 | // [[/details]] | ||
| 61 | |||
| 62 | // | ||
| 63 | // ## 1. The problem class | ||
| 64 | // | ||
| 65 | // The problem class implements the boundary conditions. It also provides | ||
| 66 | // an interface that is used by the local residual (see Part 1) to obtain the diffusion | ||
| 67 | // coefficient. The value is read from the parameter configuration tree. | ||
| 68 | // [[content]] | ||
| 69 | namespace Dumux { | ||
| 70 | template<class TypeTag> | ||
| 71 | 3 | class DiffusionTestProblem : public FVProblem<TypeTag> | |
| 72 | { | ||
| 73 | // [[details]] boilerplate code | ||
| 74 | using ParentType = FVProblem<TypeTag>; | ||
| 75 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 76 | using GlobalPosition = typename GridGeometry::LocalView::Element::Geometry::GlobalCoordinate; | ||
| 77 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 78 | using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; | ||
| 79 | using NumEqVector = Dumux::NumEqVector<PrimaryVariables>; | ||
| 80 | using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; | ||
| 81 | // [[/details]] | ||
| 82 | // In the constructor, we read the diffusion coefficient constant from the | ||
| 83 | // parameter tree (which is initialized with the content of `params.input`). | ||
| 84 | public: | ||
| 85 | 3 | DiffusionTestProblem(std::shared_ptr<const GridGeometry> gridGeometry) | |
| 86 |
6/16✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 3 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 3 times.
✓ Branch 15 taken 3 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
|
9 | : ParentType(gridGeometry) |
| 87 | { | ||
| 88 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | diffusionCoefficient_ = getParam<Scalar>("Problem.DiffusionCoefficient"); |
| 89 | 3 | } | |
| 90 | |||
| 91 | // As boundary conditions, we specify Neumann everywhere. This means, we | ||
| 92 | // have to prescribe a flux at each boundary sub control volume face | ||
| 93 | ✗ | BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const | |
| 94 | { | ||
| 95 |
1/2✓ Branch 0 taken 4488 times.
✗ Branch 1 not taken.
|
4488 | BoundaryTypes values; |
| 96 |
1/2✓ Branch 0 taken 4488 times.
✗ Branch 1 not taken.
|
4488 | values.setAllNeumann(); |
| 97 | ✗ | return values; | |
| 98 | } | ||
| 99 | |||
| 100 | // We prescribe zero flux over all of the boundary | ||
| 101 | ✗ | NumEqVector neumannAtPos(const GlobalPosition& globalPos) const | |
| 102 | 5280 | { return NumEqVector(0.0); } | |
| 103 | |||
| 104 | // The diffusion coefficient interface is used in the local residual (see Part 1). | ||
| 105 | // We can name this interface however we want as long as we adapt the calling site | ||
| 106 | // in the `LocalResidual` class in `model.hh`. | ||
| 107 | ✗ | Scalar diffusionCoefficient() const | |
| 108 | ✗ | { return diffusionCoefficient_; } | |
| 109 | |||
| 110 | private: | ||
| 111 | Scalar diffusionCoefficient_; | ||
| 112 | }; | ||
| 113 | } // end namespace Dumux | ||
| 114 | // [[/content]] | ||
| 115 | |||
| 116 | // | ||
| 117 | // ## 2. Property tag and specializations | ||
| 118 | // | ||
| 119 | // We create a new tag `DiffusionTest` that inherits all properties | ||
| 120 | // specialized for the tag `DiffusionModel` (we created this in Part 1) | ||
| 121 | // and the tag `BoxModel` which provides types relevant for the spatial | ||
| 122 | // discretization scheme (see [dumux/discretization/box.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/discretization/box.hh)). | ||
| 123 | // | ||
| 124 | // Here we choose a short form of specializing properties. The property | ||
| 125 | // system also recognizes an alias (`using`) with the property name being | ||
| 126 | // a member of the specified type tag. Note that we could also use the same mechanism | ||
| 127 | // as in (Part 1), for example: | ||
| 128 | // ```code | ||
| 129 | // template<class TypeTag> | ||
| 130 | // struct Scalar<TypeTag, TTag::DiffusionTest> | ||
| 131 | // { using type = double; }; | ||
| 132 | //``` | ||
| 133 | // which has the same effect as having an alias `Scalar = double;` | ||
| 134 | // as member of the type tag `DiffusionTest`. | ||
| 135 | // This mechanism allows for a terser code expression. | ||
| 136 | // In case both definitions are present for the same type tag, the explicit | ||
| 137 | // specialization (long form) takes precedence. | ||
| 138 | // | ||
| 139 | // [[content]] | ||
| 140 | namespace Dumux::Properties::TTag { | ||
| 141 | |||
| 142 | struct DiffusionTest | ||
| 143 | { | ||
| 144 | using InheritsFrom = std::tuple<DiffusionModel, BoxModel>; | ||
| 145 | |||
| 146 | using Scalar = double; | ||
| 147 | using Grid = Dune::YaspGrid<2>; | ||
| 148 | |||
| 149 | template<class TypeTag> | ||
| 150 | using Problem = DiffusionTestProblem<TypeTag>; | ||
| 151 | |||
| 152 | using EnableGridVolumeVariablesCache = std::true_type; | ||
| 153 | using EnableGridFluxVariablesCache = std::true_type; | ||
| 154 | using EnableGridGeometryCache = std::true_type; | ||
| 155 | }; | ||
| 156 | |||
| 157 | } // end namespace Dumux::Properties::TTag | ||
| 158 | // [[/content]] | ||
| 159 | |||
| 160 | // | ||
| 161 | // ## 3. Creating the initial solution | ||
| 162 | // | ||
| 163 | // We want to initialize the entries of the solution vector $c_{h,B}$ | ||
| 164 | // with a random field such that $c_{h,B} \sim U(0,1)$. When running | ||
| 165 | // with MPI in parallel this requires communication. With just one | ||
| 166 | // processor (`gg.gridView().comm().size() == 1`), we can just iterate | ||
| 167 | // over all entries of the solution vector and assign a uniformly | ||
| 168 | // distributed random number. However, with multiple processes, the | ||
| 169 | // solution vector only contains a subset of the degrees of freedom | ||
| 170 | // living on the processor. Moreover, some degrees of freedom exist | ||
| 171 | // on multiple processes. Each processor generates their own random | ||
| 172 | // number which means we might generate different numbers for the | ||
| 173 | // same degree of freedom. Therefore, we communicate the number. | ||
| 174 | // The idea is that each process adds its rank number (processes are | ||
| 175 | // numbered starting from 0) to its value. Then we take the minimum | ||
| 176 | // over all processes which will take return the value of the processor | ||
| 177 | // with the lowest rank. Afterwards, we subtract the added rank offset. | ||
| 178 | // | ||
| 179 | // [[content]] | ||
| 180 | // [[codeblock]] | ||
| 181 | template<class SolutionVector, class GridGeometry> | ||
| 182 | 3 | SolutionVector createInitialSolution(const GridGeometry& gg) | |
| 183 | { | ||
| 184 | 6 | SolutionVector sol(gg.numDofs()); | |
| 185 | |||
| 186 | // Generate random number and add processor offset | ||
| 187 | // For sequential run, `rank` always returns `0`. | ||
| 188 | 3 | std::mt19937 gen(0); // seed is 0 for deterministic results | |
| 189 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | Dumux::SimpleUniformDistribution<> dis(0.0, 1.0); |
| 190 | |||
| 191 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
|
6 | const auto rank = gg.gridView().comm().rank(); |
| 192 |
2/2✓ Branch 0 taken 105 times.
✓ Branch 1 taken 3 times.
|
108 | for (int n = 0; n < sol.size(); ++n) |
| 193 | 105 | sol[n] = dis(gen) + rank; | |
| 194 | |||
| 195 | // We take the value of the processor with the minimum rank | ||
| 196 | // and subtract the rank offset | ||
| 197 |
6/8✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 1 times.
|
6 | if (gg.gridView().comm().size() > 1) |
| 198 | { | ||
| 199 | Dumux::VectorCommDataHandleMin< | ||
| 200 | typename GridGeometry::VertexMapper, | ||
| 201 | SolutionVector, | ||
| 202 | GridGeometry::GridView::dimension | ||
| 203 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | > minHandle(gg.vertexMapper(), sol); |
| 204 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | gg.gridView().communicate(minHandle, Dune::All_All_Interface, Dune::ForwardCommunication); |
| 205 | |||
| 206 | // remove processor offset | ||
| 207 |
2/2✓ Branch 0 taken 56 times.
✓ Branch 1 taken 2 times.
|
58 | for (int n = 0; n < sol.size(); ++n) |
| 208 | 112 | sol[n][0] -= std::floor(sol[n][0]); | |
| 209 | } | ||
| 210 | |||
| 211 | 3 | return sol; | |
| 212 | } | ||
| 213 | // [[/codeblock]] | ||
| 214 | // [[/content]] | ||
| 215 | // | ||
| 216 | // ## 4. The main program | ||
| 217 | // | ||
| 218 | // [[content]] | ||
| 219 | 3 | int main(int argc, char** argv) | |
| 220 | { | ||
| 221 | 3 | using namespace Dumux; | |
| 222 | |||
| 223 | // First, we initialize MPI and the multithreading backend. | ||
| 224 | // This convenience function takes care that everything is setup in the right order and | ||
| 225 | // has to be called for every Dumux simulation. `Dumux::initialize` also respects | ||
| 226 | // the environment variable `DUMUX_NUM_THREADS` to restrict to amount of available cores | ||
| 227 | // for multi-threaded code parts (for example the assembly). | ||
| 228 | 3 | Dumux::initialize(argc, argv); | |
| 229 | |||
| 230 | // We initialize parameter tree including command line arguments. | ||
| 231 | // This will, per default, read all parameters from the configuration file `params.input` | ||
| 232 | // if such as file exists. Then it will look for command line arguments. For example | ||
| 233 | // `./example_diffusion -TimeLoop.TEnd 10` will set the end time to $10$ seconds. | ||
| 234 | // Command line arguments overwrite settings in the parameter file. | ||
| 235 |
9/24✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 3 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 18 taken 3 times.
✓ Branch 19 taken 3 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 3 times.
✗ 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 28 not taken.
|
12 | Parameters::init(argc, argv); |
| 236 | |||
| 237 | // We specify an alias for the model type tag. | ||
| 238 | // We will configure the assembler with this type tag that | ||
| 239 | // we specialized all these properties for above and in the model definition (Part 1). | ||
| 240 | // We can extract type information through properties specialized for the type tag | ||
| 241 | // using `GetPropType`. | ||
| 242 | 3 | using TypeTag = Properties::TTag::DiffusionTest; | |
| 243 | |||
| 244 | 3 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | |
| 245 | 3 | using Grid = GetPropType<TypeTag, Properties::Grid>; | |
| 246 | 3 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | |
| 247 | 3 | using Problem = GetPropType<TypeTag, Properties::Problem>; | |
| 248 | 3 | using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; | |
| 249 | 3 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | |
| 250 | |||
| 251 | // First, we initialize the grid. Grid parameters are taken from the input file | ||
| 252 | // from the group `[Grid]` if it exists. You can also pass any other parameter | ||
| 253 | // group (e.g. `"MyGroup"`) to `init` and then it will look in `[MyGroup.Grid]`. | ||
| 254 | 6 | GridManager<Grid> gridManager; | |
| 255 |
5/12✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
6 | gridManager.init(); |
| 256 | |||
| 257 | // We, create the finite volume grid geometry from the (leaf) grid view, | ||
| 258 | // the problem for the boundary conditions, a solution vector and a grid variables instance. | ||
| 259 |
3/8✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
6 | auto gridGeometry = std::make_shared<GridGeometry>(gridManager.grid().leafGridView()); |
| 260 |
2/6✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
6 | auto problem = std::make_shared<Problem>(gridGeometry); |
| 261 |
3/8✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
9 | auto sol = createInitialSolution<SolutionVector>(*gridGeometry); |
| 262 |
2/6✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
6 | auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry); |
| 263 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
6 | gridVariables->init(sol); |
| 264 | |||
| 265 | // We initialize the VTK output module and write out the initial concentration field | ||
| 266 |
8/16✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 3 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 3 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 3 times.
✗ Branch 22 not taken.
|
12 | VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, sol, problem->name()); |
| 267 |
7/18✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 3 times.
✓ Branch 17 taken 3 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
|
9 | vtkWriter.addVolumeVariable([](const auto& vv){ return vv.priVar(0); }, "c"); |
| 268 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | vtkWriter.write(0.0); |
| 269 | |||
| 270 | // We instantiate time loop using start and end time as well as | ||
| 271 | // the time step size from the parameter tree (`params.input`) | ||
| 272 | 3 | auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>( | |
| 273 | ✗ | getParam<Scalar>("TimeLoop.TStart", 0.0), | |
| 274 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | getParam<Scalar>("TimeLoop.Dt"), |
| 275 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
3 | getParam<Scalar>("TimeLoop.TEnd") |
| 276 |
1/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
6 | ); |
| 277 | |||
| 278 | // Next, we choose the type of assembler, linear solver and PDE solver | ||
| 279 | // and construct instances of these classes. | ||
| 280 | 3 | using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>; | |
| 281 | 3 | using LinearSolver = SSORCGIstlSolver<LinearSolverTraits<GridGeometry>, LinearAlgebraTraitsFromAssembler<Assembler>>; | |
| 282 | 3 | using Solver = Dumux::LinearPDESolver<Assembler, LinearSolver>; | |
| 283 | |||
| 284 |
2/6✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
6 | auto oldSol = sol; // copy the vector to store state of previous time step |
| 285 |
2/6✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
6 | auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, oldSol); |
| 286 |
6/14✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 3 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
|
18 | auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper()); |
| 287 |
8/20✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 3 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 3 times.
✓ Branch 19 taken 3 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 3 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
|
15 | Solver solver(assembler, linearSolver); |
| 288 | |||
| 289 | // We tell the assembler to create the system matrix and vector objects. Then we | ||
| 290 | // assemble the system matrix once and then tell the solver to reuse in every time step. | ||
| 291 | // Since we have a linear problem, the matrix is not going to change. | ||
| 292 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
6 | assembler->setLinearSystem(); |
| 293 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
6 | assembler->assembleJacobian(sol); |
| 294 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | solver.reuseMatrix(); |
| 295 | |||
| 296 | // The time loop is where most of the actual computations happen. | ||
| 297 | // We assemble and solve the linear system of equations, update the solution, | ||
| 298 | // write the solution to a VTK file and continue until the we reach the | ||
| 299 | // final simulation time. | ||
| 300 | // | ||
| 301 | // [[codeblock]] | ||
| 302 |
2/4✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
|
6 | timeLoop->start(); do |
| 303 | { | ||
| 304 | // assemble & solve | ||
| 305 |
1/2✓ Branch 1 taken 150 times.
✗ Branch 2 not taken.
|
150 | solver.solve(sol); |
| 306 | |||
| 307 | // make the new solution the old solution | ||
| 308 |
1/2✓ Branch 1 taken 150 times.
✗ Branch 2 not taken.
|
150 | oldSol = sol; |
| 309 |
2/4✓ Branch 1 taken 150 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 150 times.
✗ Branch 5 not taken.
|
300 | gridVariables->advanceTimeStep(); |
| 310 | |||
| 311 | // advance to the time loop to the next step | ||
| 312 |
2/4✓ Branch 1 taken 150 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 150 times.
✗ Branch 5 not taken.
|
300 | timeLoop->advanceTimeStep(); |
| 313 | |||
| 314 | // write VTK output (writes out the concentration field) | ||
| 315 |
3/6✓ Branch 1 taken 150 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 150 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 150 times.
✗ Branch 8 not taken.
|
450 | vtkWriter.write(timeLoop->time()); |
| 316 | |||
| 317 | // report statistics of this time step | ||
| 318 |
2/4✓ Branch 1 taken 150 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 150 times.
✗ Branch 5 not taken.
|
300 | timeLoop->reportTimeStep(); |
| 319 | |||
| 320 |
4/6✓ Branch 1 taken 150 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 150 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 147 times.
✓ Branch 7 taken 3 times.
|
300 | } while (!timeLoop->finished()); |
| 321 | |||
| 322 |
5/10✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 3 times.
✓ Branch 9 taken 3 times.
✗ Branch 10 not taken.
|
12 | timeLoop->finalize(gridGeometry->gridView().comm()); |
| 323 | |||
| 324 | 3 | return 0; | |
| 325 | } | ||
| 326 | // [[/codeblock]] | ||
| 327 | // [[/content]] | ||
| 328 |