| 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 | // # Cahn-Hilliard equation main program | ||
| 9 | // | ||
| 10 | // The file [`main.cc`](main.cc) contains the problem class `CahnHilliardTestProblem`, | ||
| 11 | // properties and traits specific to the test case as well as the `main` function. | ||
| 12 | // The setup consists of four steps: | ||
| 13 | // 1. Define the problem setting boundary conditions and the diffusion coefficient | ||
| 14 | // 2. Configure the property system reusing the model defined in Part 1 | ||
| 15 | // 3. Define a function for setting the random initial condition | ||
| 16 | // 4. The main program defining all steps of the program | ||
| 17 | // | ||
| 18 | // __Table of contents__ | ||
| 19 | // | ||
| 20 | // [TOC] | ||
| 21 | // | ||
| 22 | // We start in `main.cc` with the necessary header includes: | ||
| 23 | // [[details]] includes | ||
| 24 | #include <config.h> | ||
| 25 | #include <type_traits> | ||
| 26 | #include <random> | ||
| 27 | |||
| 28 | #include <dumux/common/initialize.hh> | ||
| 29 | #include <dumux/common/properties.hh> | ||
| 30 | #include <dumux/common/parameters.hh> | ||
| 31 | #include <dumux/common/numeqvector.hh> | ||
| 32 | #include <dumux/common/boundarytypes.hh> | ||
| 33 | #include <dumux/common/fvproblem.hh> | ||
| 34 | |||
| 35 | #include <dumux/discretization/box.hh> | ||
| 36 | #include <dumux/linear/linearsolvertraits.hh> | ||
| 37 | #include <dumux/linear/linearalgebratraits.hh> | ||
| 38 | #include <dumux/linear/istlsolvers.hh> | ||
| 39 | #include <dumux/nonlinear/newtonsolver.hh> | ||
| 40 | #include <dumux/assembly/fvassembler.hh> | ||
| 41 | |||
| 42 | #include <dune/grid/yaspgrid.hh> | ||
| 43 | #include <dumux/io/chrono.hh> | ||
| 44 | #include <dumux/io/vtkoutputmodule.hh> | ||
| 45 | #include <dumux/io/grid/gridmanager_yasp.hh> | ||
| 46 | |||
| 47 | #include "model.hh" | ||
| 48 | // [[/details]] | ||
| 49 | |||
| 50 | // | ||
| 51 | // ## 1. The problem class | ||
| 52 | // | ||
| 53 | // The problem class implements boundary conditions and the source term. | ||
| 54 | // It also provides an interface that is used by the local residual (see Part 1) to obtain | ||
| 55 | // the mobility and surface tension. The values are read from the parameter configuration tree. | ||
| 56 | // [[content]] | ||
| 57 | namespace Dumux { | ||
| 58 | template<class TypeTag> | ||
| 59 | 3 | class CahnHilliardTestProblem : public FVProblem<TypeTag> | |
| 60 | { | ||
| 61 | // [[details]] boilerplate code | ||
| 62 | using ParentType = FVProblem<TypeTag>; | ||
| 63 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 64 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 65 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 66 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
| 67 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 68 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 69 | |||
| 70 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 71 | using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; | ||
| 72 | using NumEqVector = Dumux::NumEqVector<PrimaryVariables>; | ||
| 73 | using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; | ||
| 74 | using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; | ||
| 75 | // [[/details]] | ||
| 76 | // In the constructor, we read the parameter constants from the | ||
| 77 | // parameter tree (which is initialized with the content of `params.input`). | ||
| 78 | public: | ||
| 79 | 3 | CahnHilliardTestProblem(std::shared_ptr<const GridGeometry> gridGeometry) | |
| 80 |
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) |
| 81 | { | ||
| 82 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | mobility_ = getParam<Scalar>("Problem.Mobility"); |
| 83 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | surfaceTension_ = getParam<Scalar>("Problem.SurfaceTension"); |
| 84 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | energyScale_ = getParam<Scalar>("Problem.EnergyScale"); |
| 85 | 3 | } | |
| 86 | |||
| 87 | // In the `source` function, we implement the derivative of the free energy. | ||
| 88 | // This demonstrates how parts of the local residual can be split into model specific | ||
| 89 | // parts (see `CahnHilliardModelLocalResidual`) and parts that might change from scenario to scenario. | ||
| 90 | template<class ElementVolumeVariables> | ||
| 91 | 3672800 | NumEqVector source(const Element &element, | |
| 92 | const FVElementGeometry& fvGeometry, | ||
| 93 | const ElementVolumeVariables& elemVolVars, | ||
| 94 | const SubControlVolume &scv) const | ||
| 95 | { | ||
| 96 | 3672800 | NumEqVector values(0.0); | |
| 97 | 11018400 | const auto& c = elemVolVars[scv].concentration(); | |
| 98 | 7345600 | values[Indices::chemicalPotentialEqIdx] = -energyScale_*2.0*c*(2.0*c*c - 3*c + 1); | |
| 99 | 3672800 | return values; | |
| 100 | } | ||
| 101 | |||
| 102 | // We choose boundary flux (or Neumann) conditions for all equations on the entire boundary, | ||
| 103 | // while specifying zero flux for both equations. | ||
| 104 | // [[codeblock]] | ||
| 105 | ✗ | BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const | |
| 106 | { | ||
| 107 | 72560 | BoundaryTypes values; | |
| 108 |
1/2✓ Branch 0 taken 72560 times.
✗ Branch 1 not taken.
|
72560 | values.setAllNeumann(); |
| 109 | ✗ | return values; | |
| 110 | } | ||
| 111 | |||
| 112 | ✗ | NumEqVector neumannAtPos(const GlobalPosition& globalPos) const | |
| 113 | 1335104 | { return { 0.0, 0.0 }; } | |
| 114 | // [[/codeblock]] | ||
| 115 | |||
| 116 | // The parameters interfaces are used in the local residual (see Part 1). | ||
| 117 | // We can name this interface however we want as long as we adapt the calling site | ||
| 118 | // in the `CahnHilliardModelLocalResidual` class in `model.hh`. | ||
| 119 | // [[codeblock]] | ||
| 120 | ✗ | Scalar mobility() const | |
| 121 | ✗ | { return mobility_; } | |
| 122 | |||
| 123 | ✗ | Scalar surfaceTension() const | |
| 124 | ✗ | { return surfaceTension_; } | |
| 125 | |||
| 126 | private: | ||
| 127 | Scalar mobility_; | ||
| 128 | Scalar surfaceTension_; | ||
| 129 | Scalar energyScale_; | ||
| 130 | }; | ||
| 131 | } // end namespace Dumux | ||
| 132 | // [[/codeblock]] | ||
| 133 | // [[/content]] | ||
| 134 | |||
| 135 | // | ||
| 136 | // ## 2. Property tag and specializations | ||
| 137 | // | ||
| 138 | // We create a new tag `DiffusionTest` that inherits all properties | ||
| 139 | // specialized for the tag `DiffusionModel` (we created this in Part 1) | ||
| 140 | // and the tag `BoxModel` which provides types relevant for the spatial | ||
| 141 | // discretization scheme (see [dumux/discretization/box.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/discretization/box.hh)). | ||
| 142 | // | ||
| 143 | // [[content]] | ||
| 144 | namespace Dumux::Properties::TTag { | ||
| 145 | struct CahnHilliardTest | ||
| 146 | { | ||
| 147 | using InheritsFrom = std::tuple<CahnHilliardModel, BoxModel>; | ||
| 148 | |||
| 149 | using Grid = Dune::YaspGrid<2>; | ||
| 150 | |||
| 151 | template<class TypeTag> | ||
| 152 | using Problem = CahnHilliardTestProblem<TypeTag>; | ||
| 153 | |||
| 154 | using EnableGridGeometryCache = std::true_type; | ||
| 155 | using EnableGridVolumeVariablesCache = std::true_type; | ||
| 156 | using EnableGridFluxVariablesCache = std::true_type; | ||
| 157 | }; | ||
| 158 | } // end namespace Dumux::Properties | ||
| 159 | // [[/content]] | ||
| 160 | // | ||
| 161 | // ## 3. Creating the initial solution | ||
| 162 | // | ||
| 163 | // To initialize with a random field in parallel, where each processor | ||
| 164 | // creates its own random number sequence, we need to communicate the | ||
| 165 | // resulting values on the processor border and overlap. | ||
| 166 | // See [Diffusion equation example](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/examples/diffusion/docs/main.md). | ||
| 167 | // for details. Here in addition, we need to provide a custom scatter operation | ||
| 168 | // that handles vector types. We only need to communicate the first entry (concentration). | ||
| 169 | // | ||
| 170 | // [[content]] | ||
| 171 | // [[codeblock]] | ||
| 172 | // functor for data communication with MPI | ||
| 173 | struct MinScatter | ||
| 174 | { | ||
| 175 | template<class A, class B> | ||
| 176 | static void apply(A& a, const B& b) | ||
| 177 | ✗ | { a[0] = std::min(a[0], b[0]); } | |
| 178 | }; | ||
| 179 | |||
| 180 | // create the random initial solution | ||
| 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 runs `rank` always returns `0`. | ||
| 188 | 3 | std::mt19937 gen(0); // seed is 0 for deterministic results | |
| 189 | 3 | std::uniform_real_distribution<> dis(0.0, 1.0); | |
| 190 |
2/2✓ Branch 0 taken 297 times.
✓ Branch 1 taken 3 times.
|
300 | for (int n = 0; n < sol.size(); ++n) |
| 191 | { | ||
| 192 |
2/4✗ Branch 1 not taken.
✓ Branch 2 taken 297 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 297 times.
|
297 | sol[n][0] = 0.42 + 0.02*(0.5-dis(gen)) + gg.gridView().comm().rank(); |
| 193 | 594 | sol[n][1] = 0.0; | |
| 194 | } | ||
| 195 | |||
| 196 | // We take the value of the processor with the minimum rank | ||
| 197 | // and subtract the rank offset | ||
| 198 |
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) |
| 199 | { | ||
| 200 | Dumux::VectorCommDataHandle< | ||
| 201 | typename GridGeometry::VertexMapper, | ||
| 202 | SolutionVector, | ||
| 203 | GridGeometry::GridView::dimension, | ||
| 204 | MinScatter | ||
| 205 |
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); |
| 206 |
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); |
| 207 | |||
| 208 | // Remove processor offset | ||
| 209 |
2/2✓ Branch 0 taken 176 times.
✓ Branch 1 taken 2 times.
|
178 | for (int n = 0; n < sol.size(); ++n) |
| 210 | 704 | sol[n][0] -= std::floor(sol[n][0]); | |
| 211 | } | ||
| 212 | 3 | return sol; | |
| 213 | } | ||
| 214 | // [[/codeblock]] | ||
| 215 | // [[/content]] | ||
| 216 | // | ||
| 217 | // ## 4. The main program | ||
| 218 | // | ||
| 219 | // The `main` function sets up the simulation framework, initializes runtime parameters, | ||
| 220 | // creates the grid and storage vectors | ||
| 221 | // for the variables, primary and secondary. It specifies and constructs and assembler, which | ||
| 222 | // assembles the discretized residual and system matrix (Jacobian of the model residual), as well as | ||
| 223 | // a linear solver. A Newton method is used to solve the nonlinear equations where in each Newton iteration | ||
| 224 | // the Jacobian and the residual needs to be reassembled and the resulting linear system is solved. | ||
| 225 | // The time loop controls the time stepping. Here, we let the Newton solver suggest an adaption of | ||
| 226 | // the time step size based on the convergence history of the nonlinear solver. | ||
| 227 | // | ||
| 228 | // [[content]] | ||
| 229 | 3 | int main(int argc, char** argv) | |
| 230 | { | ||
| 231 | 3 | using namespace Dumux; | |
| 232 | |||
| 233 | // We initialize MPI and/or multithreading backend. When not running | ||
| 234 | // in parallel the MPI setup is skipped. | ||
| 235 | 3 | Dumux::initialize(argc, argv); | |
| 236 | |||
| 237 | // Then we initialize parameter tree. | ||
| 238 |
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); |
| 239 | |||
| 240 | // We create an alias for the type tag for this problem | ||
| 241 | // and extract type information through properties. | ||
| 242 | 3 | using TypeTag = Properties::TTag::CahnHilliardTest; | |
| 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 | // We initialize the grid, grid geometry, problem, solution vector, and grid variables. | ||
| 252 | 6 | GridManager<Grid> gridManager; | |
| 253 |
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(); |
| 254 | |||
| 255 |
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()); |
| 256 |
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); |
| 257 |
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); |
| 258 |
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); |
| 259 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
6 | gridVariables->init(sol); |
| 260 | |||
| 261 | // We initialize the VTK output module and write out the initial concentration and chemical potential | ||
| 262 |
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()); |
| 263 |
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.concentration(); }, "c"); |
| 264 |
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.chemicalPotential(); }, "mu"); |
| 265 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | vtkWriter.write(0.0); |
| 266 | |||
| 267 | // We instantiate time loop using start and end time as well as | ||
| 268 | // the time step size from the parameter tree (`params.input`) | ||
| 269 | 3 | auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>( | |
| 270 |
1/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
3 | Chrono::toSeconds(getParam("TimeLoop.TStart", "0")), |
| 271 |
3/6✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 3 times.
|
6 | Chrono::toSeconds(getParam("TimeLoop.InitialTimeStepSize")), |
| 272 |
4/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 9 not taken.
✓ Branch 10 taken 3 times.
|
6 | Chrono::toSeconds(getParam("TimeLoop.TEnd")) |
| 273 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
|
6 | ); |
| 274 | |||
| 275 | // We set the maximum time step size allowed in the adaptive time stepping scheme. | ||
| 276 |
5/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 12 not taken.
✓ Branch 13 taken 3 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
|
6 | timeLoop->setMaxTimeStepSize(Chrono::toSeconds(getParam("TimeLoop.MaxTimeStepSize"))); |
| 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 = SSORBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>, LinearAlgebraTraitsFromAssembler<Assembler>>; | |
| 282 | 3 | using Solver = Dumux::NewtonSolver<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; |
| 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 18 taken 3 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 3 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
|
15 | Solver solver(assembler, linearSolver); |
| 288 | |||
| 289 | // The time loop is where most of the actual computations happen. | ||
| 290 | // We assemble and solve the linear system of equations, update the solution, | ||
| 291 | // write the solution to a VTK file and continue until the we reach the | ||
| 292 | // final simulation time. | ||
| 293 | // | ||
| 294 | // [[codeblock]] | ||
| 295 |
2/4✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
|
6 | timeLoop->start(); do |
| 296 | { | ||
| 297 | // Assemble & solve | ||
| 298 | // Passing the time loop enables the solver to repeat the time step | ||
| 299 | // with a reduced time step size if the Newton solver does not converge. | ||
| 300 |
2/4✓ Branch 1 taken 306 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 306 times.
✗ Branch 5 not taken.
|
612 | solver.solve(sol, *timeLoop); |
| 301 | |||
| 302 | // make the new solution the old solution | ||
| 303 |
1/2✓ Branch 1 taken 306 times.
✗ Branch 2 not taken.
|
306 | oldSol = sol; |
| 304 |
2/4✓ Branch 1 taken 306 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 306 times.
✗ Branch 5 not taken.
|
612 | gridVariables->advanceTimeStep(); |
| 305 | |||
| 306 | // advance to the time loop to the next step | ||
| 307 |
2/4✓ Branch 1 taken 306 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 306 times.
✗ Branch 5 not taken.
|
612 | timeLoop->advanceTimeStep(); |
| 308 | |||
| 309 | // write VTK output (concentration field and chemical potential) | ||
| 310 |
3/6✓ Branch 1 taken 306 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 306 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 306 times.
✗ Branch 8 not taken.
|
918 | vtkWriter.write(timeLoop->time()); |
| 311 | |||
| 312 | // report statistics of this time step | ||
| 313 |
2/4✓ Branch 1 taken 306 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 306 times.
✗ Branch 5 not taken.
|
612 | timeLoop->reportTimeStep(); |
| 314 | |||
| 315 | // set new dt as suggested by the Newton solver | ||
| 316 |
5/10✗ Branch 0 not taken.
✓ Branch 1 taken 306 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 306 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 306 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 306 times.
✓ Branch 9 taken 306 times.
✗ Branch 10 not taken.
|
1530 | timeLoop->setTimeStepSize(solver.suggestTimeStepSize(timeLoop->timeStepSize())); |
| 317 | |||
| 318 |
4/6✓ Branch 1 taken 306 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 306 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 303 times.
✓ Branch 7 taken 3 times.
|
612 | } while (!timeLoop->finished()); |
| 319 | |||
| 320 | // print the final report | ||
| 321 |
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()); |
| 322 | 3 | return 0; | |
| 323 | } | ||
| 324 | // [[/codeblock]] | ||
| 325 | // [[/content]] | ||
| 326 |