| 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 | // # The main file (`main.cc`) | ||
| 8 | // | ||
| 9 | // This file contains the `main` function and implements the main program | ||
| 10 | // flow. We initialize the simulation framework, initialize parameters, | ||
| 11 | // create the grids (using parameters from the configuration file `params.input`), | ||
| 12 | // create vector to store the primary and secondary variables, construct an | ||
| 13 | // assembler that can assembly the the residual (discrete equations) and | ||
| 14 | // the system matrix (Jacobian of the residual), and create a linear solver | ||
| 15 | // that solves the linear system arising in each time step. The time loop | ||
| 16 | // implements a simple time stepping scheme with constant time step size. | ||
| 17 | // | ||
| 18 | // [[content]] | ||
| 19 | // | ||
| 20 | // ### Included header files | ||
| 21 | // [[details]] includes | ||
| 22 | // [[codeblock]] | ||
| 23 | #include <config.h> | ||
| 24 | |||
| 25 | #include <vector> | ||
| 26 | #include <numeric> | ||
| 27 | |||
| 28 | #include <dune/common/timer.hh> | ||
| 29 | #include <dune/common/fvector.hh> | ||
| 30 | |||
| 31 | #include <dumux/common/properties.hh> | ||
| 32 | #include <dumux/common/parameters.hh> | ||
| 33 | #include <dumux/common/initialize.hh> | ||
| 34 | |||
| 35 | #include <dumux/io/vtkoutputmodule.hh> | ||
| 36 | #include <dumux/io/container.hh> | ||
| 37 | #include <dumux/io/format.hh> | ||
| 38 | #include <dumux/io/grid/gridmanager_yasp.hh> | ||
| 39 | #include <dumux/io/grid/gridmanager_foam.hh> | ||
| 40 | |||
| 41 | #include <dumux/multidomain/traits.hh> | ||
| 42 | #include <dumux/multidomain/fvassembler.hh> | ||
| 43 | #include <dumux/multidomain/newtonsolver.hh> | ||
| 44 | |||
| 45 | #include "bloodflow.hh" | ||
| 46 | #include "solver.hh" | ||
| 47 | #include "properties.hh" | ||
| 48 | // [[/codeblock]] | ||
| 49 | // [[/details]] | ||
| 50 | // | ||
| 51 | // ## The main function | ||
| 52 | // We will now discuss the main program flow implemented within the `main` function. | ||
| 53 | // At the beginning of each program using Dune, an instance of `Dune::MPIHelper` has to | ||
| 54 | // be created. Moreover, we parse the run-time arguments from the command line and the | ||
| 55 | // input file: | ||
| 56 | 1 | int main(int argc, char** argv) { | |
| 57 | // enable all symbols of the Dumux namespace without having to prefix with ``Dumux::`` | ||
| 58 | 1 | using namespace Dumux; | |
| 59 | // | ||
| 60 | // ### initialization and parameter setup | ||
| 61 | // [[codeblock]] | ||
| 62 | // initialize (e.g. multi-threading backend) | ||
| 63 | 1 | initialize(argc, argv); | |
| 64 | |||
| 65 | // parse command line arguments and input file | ||
| 66 |
9/24✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 1 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.
|
6 | Parameters::init(argc, argv); |
| 67 | // [[/codeblock]] | ||
| 68 | |||
| 69 | // Define the sub problem type tags (collection of model properties) | ||
| 70 | 1 | using Model3D = Properties::TTag::TissueTransportModel; | |
| 71 | 1 | using Model1D = Properties::TTag::NetworkTransportModel; | |
| 72 | |||
| 73 | // Start the timer (to measure code execution time) | ||
| 74 | 1 | Dune::Timer timer; | |
| 75 | |||
| 76 | // ### Grid setup | ||
| 77 | // [[codeblock]] | ||
| 78 | // Create grids (see input file groups "[Tissue.Grid]" and "[Network.Grid]") | ||
| 79 | // See properties.hh for the types of Grid | ||
| 80 | 1 | using Grid3D = GetPropType<Model3D, Properties::Grid>; | |
| 81 | 2 | GridManager<Grid3D> gridManager3D; | |
| 82 |
5/12✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
2 | gridManager3D.init("Tissue"); |
| 83 | |||
| 84 | 1 | using Grid1D = GetPropType<Model1D, Properties::Grid>; | |
| 85 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
2 | GridManager<Grid1D> gridManager1D; |
| 86 |
5/12✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
2 | gridManager1D.init("Network"); |
| 87 | |||
| 88 | // get data from grid file like vessel radii and fluxes in the network | ||
| 89 |
1/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
2 | auto gridData1D = gridManager1D.getGridData(); |
| 90 | |||
| 91 | // we compute on the leaf grid view (the refined elements) | ||
| 92 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
1 | const auto& gridView3D = gridManager3D.grid().leafGridView(); |
| 93 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
1 | const auto& gridView1D = gridManager1D.grid().leafGridView(); |
| 94 | |||
| 95 | // create the finite volume grid geometries on top of the grid | ||
| 96 | // this enables us to iterate over (sub-)control-volumes and faces | ||
| 97 | 1 | using GridGeometry3D = GetPropType<Model3D, Properties::GridGeometry>; | |
| 98 | 1 | using GridGeometry1D = GetPropType<Model1D, Properties::GridGeometry>; | |
| 99 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
2 | auto gridGeometry3D = std::make_shared<GridGeometry3D>(gridView3D); |
| 100 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
2 | auto gridGeometry1D = std::make_shared<GridGeometry1D>(gridView1D); |
| 101 | |||
| 102 | // print info about the grid setup time | ||
| 103 | 1 | const auto gridSetupTime = timer.elapsed(); | |
| 104 |
2/4✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
|
3 | std::cout << "Setting up grid geometry took " << gridSetupTime << " seconds." << std::endl; |
| 105 | // [[/codeblock]] | ||
| 106 | // | ||
| 107 | // ### Coupling manager | ||
| 108 | // | ||
| 109 | // [[codeblock]] | ||
| 110 | // the mixed dimension type traits (create coupled model properties) | ||
| 111 | 1 | using CoupledModelTraits = MultiDomainTraits<Model3D, Model1D>; | |
| 112 | 1 | constexpr auto domain3D = CoupledModelTraits::template SubDomain<0>::Index(); | |
| 113 | 1 | constexpr auto domain1D = CoupledModelTraits::template SubDomain<1>::Index(); | |
| 114 | |||
| 115 | // the coupling manager (handles data mapping between 1D and 3D domains) | ||
| 116 | 1 | using CouplingManager = Properties::CouplingTransport; | |
| 117 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
2 | auto couplingManager = std::make_shared<CouplingManager>(gridGeometry3D, gridGeometry1D); |
| 118 | |||
| 119 | // the problem (initial and boundary conditions, source and coupling terms) | ||
| 120 | // these class are defined in `problem.hh` and then defines as "properties" in `properties.hh` | ||
| 121 | // the constructor is whatever is defined as a constructor in `problem.hh`, so the user is | ||
| 122 | // free to change this to whatever works. Here we added the argument `gridData1D` to hand | ||
| 123 | // the grid parameters read from the grid file to the problem class instance | ||
| 124 | 1 | using TransportProblem3D = GetPropType<Model3D, Properties::Problem>; | |
| 125 | 1 | using TransportProblem1D = GetPropType<Model1D, Properties::Problem>; | |
| 126 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
2 | auto transportProblem3D = std::make_shared<TransportProblem3D>(gridGeometry3D, couplingManager, "Tissue"); |
| 127 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
2 | auto transportProblem1D = std::make_shared<TransportProblem1D>(gridGeometry1D, couplingManager, gridData1D, "Network"); |
| 128 | |||
| 129 | // solve blood flow problem to obtain volume fluxes | ||
| 130 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
2 | auto bloodVolumeFluxes = computeBloodVolumeFluxes(gridGeometry1D, gridData1D); |
| 131 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
3 | transportProblem1D->spatialParams().setVolumeFluxes(bloodVolumeFluxes); |
| 132 | // [[/codeblock]] | ||
| 133 | // | ||
| 134 | // ### Solution vector | ||
| 135 | // | ||
| 136 | // [[codeblock]] | ||
| 137 | // the solution vector (here: mole fractions of the transported solute) | ||
| 138 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
2 | CoupledModelTraits::SolutionVector sol; |
| 139 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
3 | transportProblem3D->applyInitialSolution(sol[domain3D]); |
| 140 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
3 | transportProblem1D->applyInitialSolution(sol[domain1D]); |
| 141 | |||
| 142 | // this is a time-dependent problem, so we create a second solution vector that | ||
| 143 | // represents the solution at the previous time level | ||
| 144 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
2 | auto oldSol = sol; |
| 145 | // [[/codeblock]] | ||
| 146 | // | ||
| 147 | // ### Coupling manager | ||
| 148 | // | ||
| 149 | // [[codeblock]] | ||
| 150 | // initialize coupling | ||
| 151 |
3/10✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
3 | couplingManager->init(transportProblem3D, transportProblem1D, sol); |
| 152 | // in the case of 1D-3D the coupling needs point sources to be enabled for both sub-problems | ||
| 153 | // each point sources constitutes an integration point for the coupling operator | ||
| 154 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | transportProblem3D->computePointSourceMap(); |
| 155 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | transportProblem1D->computePointSourceMap(); |
| 156 | |||
| 157 | // the grid variables (secondary variables and caching of computed variables) | ||
| 158 | // for instance the effective diffusion coefficient that is computed from | ||
| 159 | // the free diffusion coefficient (fluid system) and the porosity and tortuosity (spatial parameter) | ||
| 160 | 1 | using GridVariables3D = GetPropType<Model3D, Properties::GridVariables>; | |
| 161 |
1/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
2 | auto gridVariables3D = std::make_shared<GridVariables3D>(transportProblem3D, gridGeometry3D); |
| 162 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
3 | gridVariables3D->init(sol[domain3D]); |
| 163 | 1 | using GridVariables1D = GetPropType<Model1D, Properties::GridVariables>; | |
| 164 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
2 | auto gridVariables1D = std::make_shared<GridVariables1D>(transportProblem1D, gridGeometry1D); |
| 165 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
3 | gridVariables1D->init(sol[domain1D]); |
| 166 | // [[/codeblock]] | ||
| 167 | // | ||
| 168 | // ### VTK output | ||
| 169 | // | ||
| 170 | // [[codeblock]] | ||
| 171 | // initialize the VTK output | ||
| 172 | 1 | using Solution3D = std::decay_t<decltype(sol[domain3D])>; | |
| 173 |
9/18✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
|
4 | VtkOutputModule<GridVariables3D, Solution3D> vtkWriter3D(*gridVariables3D, sol[domain3D], transportProblem3D->name()); |
| 174 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | GetPropType<Model3D, Properties::IOFields>::initOutputModule(vtkWriter3D); |
| 175 | 1 | using Solution1D = std::decay_t<decltype(sol[domain1D])>; | |
| 176 |
8/16✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
|
3 | VtkOutputModule<GridVariables1D, Solution1D> vtkWriter1D(*gridVariables1D, sol[domain1D], transportProblem1D->name()); |
| 177 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | GetPropType<Model1D, Properties::IOFields>::initOutputModule(vtkWriter1D); |
| 178 | |||
| 179 | // we compute with mole fraction as primary variables because they are continuous across material interfaces | ||
| 180 | // for the output, we compute fluid concentration c = rho*x and total/voxel concentration ct = porosity*c | ||
| 181 | // mmol/l is equivalent to mol/m^3, the balance equation is in mol/s (change of tracer amount per time) | ||
| 182 |
7/18✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
|
4 | vtkWriter3D.addVolumeVariable( [](const auto& v){ return v.molarDensity()*v.moleFraction(0,0); }, "c fluid (mmol/l)"); |
| 183 |
7/18✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
|
4 | vtkWriter3D.addVolumeVariable( [](const auto& v){ return v.molarDensity()*v.moleFraction(0,0)*v.porosity(); }, "c total (mmol/l)"); |
| 184 | |||
| 185 | // for the network domain we want to output the vessel lumen radius for visualization | ||
| 186 | // and the outer vessel radius (where we couple with the tissue domain) | ||
| 187 |
8/18✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
3 | vtkWriter1D.addField(transportProblem1D->spatialParams().getVesselRadii(), "vessel radius (m)"); |
| 188 |
8/18✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
3 | vtkWriter1D.addField(transportProblem1D->spatialParams().getOuterRadii(), "outer radius (m)"); |
| 189 |
8/18✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 1 times.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
2 | vtkWriter1D.addField(transportProblem1D->spatialParams().getVesselSegmentLengths(), "length"); |
| 190 |
7/18✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
|
4 | vtkWriter1D.addVolumeVariable( [](const auto& v){ return v.molarDensity()*v.moleFraction(0,0); }, "c fluid (mmol/l)"); |
| 191 | |||
| 192 | // write out initial solution | ||
| 193 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | vtkWriter3D.write(0.0, Dune::VTK::OutputType::appendedraw); |
| 194 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | vtkWriter1D.write(0.0, Dune::VTK::OutputType::appendedraw); |
| 195 | // [[/codeblock]] | ||
| 196 | // | ||
| 197 | // ### Time loop configuration, assembler and linear solver | ||
| 198 | // | ||
| 199 | // [[codeblock]] | ||
| 200 | //////////////////////////////////////////////////////////// | ||
| 201 | // Time loop with assembler and linear solver | ||
| 202 | //////////////////////////////////////////////////////////// | ||
| 203 | |||
| 204 | // get some time loop parameters and make time loop | ||
| 205 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | const auto dt = getParam<double>("TimeLoop.DtInitial"); |
| 206 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | const auto tEnd = getParam<double>("TimeLoop.TEnd"); |
| 207 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
2 | auto timeLoop = std::make_shared<CheckPointTimeLoop<double>>(0.0, dt, tEnd); |
| 208 | |||
| 209 | // the assembler with time loop for a coupled transient problem | ||
| 210 | 1 | using Assembler = MultiDomainFVAssembler<CoupledModelTraits, CouplingManager, DiffMethod::numeric>; | |
| 211 | 1 | auto assembler = std::make_shared<Assembler>( | |
| 212 | 1 | std::make_tuple(transportProblem3D, transportProblem1D), | |
| 213 | 2 | std::make_tuple(gridGeometry3D, gridGeometry1D), | |
| 214 |
0/2✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
2 | std::make_tuple(gridVariables3D, gridVariables1D), |
| 215 | couplingManager, timeLoop, oldSol | ||
| 216 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
2 | ); |
| 217 | |||
| 218 | // construct the linear solver (see solver.hh) | ||
| 219 | 1 | using LinearSolver = Example::BlockDiagILU0BiCGSTABSolver<typename Assembler::JacobianMatrix, typename Assembler::ResidualType>; | |
| 220 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
2 | auto linearSolver = std::make_shared<LinearSolver>(); |
| 221 | // [[/codeblock]] | ||
| 222 | // | ||
| 223 | // ### Output of the tracer concentration | ||
| 224 | // | ||
| 225 | // [[codeblock]] | ||
| 226 | // output time and tracer amount in tissue | ||
| 227 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
2 | const auto filePrefix = getParam<std::string>("Problem.Name"); |
| 228 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
|
2 | std::vector<Dune::FieldVector<double, 2>> tracerAmounts; |
| 229 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | tracerAmounts.reserve(std::size_t(std::ceil(tEnd/dt))); |
| 230 |
4/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
|
4 | const auto initTracerAmount = transportProblem3D->computeTracerAmount(sol[domain3D], *gridVariables3D); |
| 231 |
3/8✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
1 | std::cout << Fmt::format("Initial amount of tracer in tissue: {:.2e}\n", initTracerAmount); |
| 232 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | tracerAmounts.push_back(Dune::FieldVector<double, 2>{ 0.0, initTracerAmount }); |
| 233 | // [[/codeblock]] | ||
| 234 | // | ||
| 235 | // ### Time loop | ||
| 236 | // | ||
| 237 | // [[codeblock]] | ||
| 238 | // start the time loop | ||
| 239 |
2/4✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
2 | timeLoop->start(); |
| 240 |
4/6✓ Branch 1 taken 101 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 101 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 100 times.
✓ Branch 7 taken 1 times.
|
202 | while (!timeLoop->finished()) |
| 241 | { | ||
| 242 | // construct some timers to measure the execution time of the time integration steps | ||
| 243 | 100 | Dune::Timer assembleTimer(false), solveTimer(false), updateTimer(false), outputTimer(false); | |
| 244 |
1/2✓ Branch 2 taken 100 times.
✗ Branch 3 not taken.
|
200 | std::cout << "\nAssemble linear system " << std::endl; |
| 245 | |||
| 246 | // assemble stiffness matrix and/or residual | ||
| 247 | // we only assemble the Jacobian (stiffness matrix for a linear problem) once | ||
| 248 | // in the subsequent time steps only the right-hand side (residual) changes | ||
| 249 |
1/2✓ Branch 0 taken 100 times.
✗ Branch 1 not taken.
|
100 | assembleTimer.start(); |
| 250 |
2/4✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
|
200 | couplingManager->updateSolution(sol); |
| 251 |
4/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 99 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 99 times.
|
200 | if (timeLoop->timeStepIndex() == 0) |
| 252 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | assembler->assembleJacobianAndResidual(sol); |
| 253 | else | ||
| 254 |
3/6✓ Branch 1 taken 99 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 99 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 99 times.
✗ Branch 8 not taken.
|
297 | assembler->assembleResidual(sol); |
| 255 | 100 | assembleTimer.stop(); | |
| 256 | |||
| 257 |
5/14✓ Branch 2 taken 100 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 100 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 100 times.
✗ Branch 9 not taken.
✓ Branch 12 taken 100 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 100 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
400 | std::cout << "Solve linear system (" << linearSolver->name() << ")" << std::endl; |
| 258 | |||
| 259 | // solve linear system | ||
| 260 | // the first time we tell the solver about the Jacobian | ||
| 261 | // then we solve with updated right-hand sides (residual) | ||
| 262 |
1/2✓ Branch 0 taken 100 times.
✗ Branch 1 not taken.
|
100 | solveTimer.start(); |
| 263 |
1/2✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
|
200 | auto deltaSol = sol; |
| 264 |
4/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 99 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 99 times.
|
200 | if (timeLoop->timeStepIndex() == 0) |
| 265 |
4/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
|
4 | linearSolver->setup(assembler->jacobian()); |
| 266 | |||
| 267 |
4/8✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 100 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 100 times.
✗ Branch 11 not taken.
|
400 | const bool converged = linearSolver->solve(deltaSol, assembler->residual()); |
| 268 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 100 times.
|
100 | if (!converged) |
| 269 | ✗ | DUNE_THROW(Dune::MathError, "Linear solver did not converge!"); | |
| 270 | 100 | solveTimer.stop(); | |
| 271 | |||
| 272 | // update solution and primary and secondary variables | ||
| 273 | // and tell the coupling manager about the new solution | ||
| 274 |
1/2✓ Branch 0 taken 100 times.
✗ Branch 1 not taken.
|
100 | updateTimer.start(); |
| 275 | 100 | sol -= deltaSol; | |
| 276 |
2/4✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
|
200 | couplingManager->updateSolution(sol); |
| 277 |
3/6✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 100 times.
✗ Branch 8 not taken.
|
300 | gridVariables3D->update(sol[domain3D]); |
| 278 |
3/6✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 100 times.
✗ Branch 8 not taken.
|
300 | gridVariables1D->update(sol[domain1D]); |
| 279 | |||
| 280 | // make the new solution the old solution | ||
| 281 |
1/2✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
|
100 | oldSol = sol; |
| 282 |
2/4✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
|
200 | gridVariables3D->advanceTimeStep(); |
| 283 |
2/4✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
|
200 | gridVariables1D->advanceTimeStep(); |
| 284 | 100 | updateTimer.stop(); | |
| 285 | |||
| 286 | // advance to the time loop to the next step | ||
| 287 |
2/4✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
|
200 | timeLoop->advanceTimeStep(); |
| 288 | |||
| 289 | // write vtk output | ||
| 290 |
1/2✓ Branch 0 taken 100 times.
✗ Branch 1 not taken.
|
100 | outputTimer.start(); |
| 291 |
3/6✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 100 times.
✗ Branch 8 not taken.
|
300 | vtkWriter3D.write(timeLoop->time(), Dune::VTK::OutputType::appendedraw); |
| 292 |
3/6✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 100 times.
✗ Branch 8 not taken.
|
300 | vtkWriter1D.write(timeLoop->time(), Dune::VTK::OutputType::appendedraw); |
| 293 | |||
| 294 | // Compute total tracer mass in domain and write to file | ||
| 295 |
4/8✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 100 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 100 times.
✗ Branch 11 not taken.
|
400 | const auto tracerAmount = transportProblem3D->computeTracerAmount(sol[domain3D], *gridVariables3D); |
| 296 |
3/8✓ Branch 2 taken 100 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 100 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 100 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
100 | std::cout << Fmt::format("Current amount of tracer in ECS: {:.2e}\n", tracerAmount); |
| 297 |
4/8✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 100 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 100 times.
✗ Branch 11 not taken.
|
400 | tracerAmounts.push_back(Dune::FieldVector<double, 2>{ timeLoop->time(), tracerAmount }); |
| 298 |
2/6✓ Branch 3 taken 100 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 100 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
200 | Dumux::writeContainerToFile(tracerAmounts, Fmt::format("{}_tracer_amounts.dat", filePrefix)); |
| 299 | 100 | outputTimer.stop(); | |
| 300 | |||
| 301 | // output timings | ||
| 302 | 100 | const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed() + outputTimer.elapsed(); | |
| 303 | 100 | std::cout << "Assemble/solve/update/output time: " | |
| 304 |
2/4✓ Branch 2 taken 100 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 100 times.
✗ Branch 8 not taken.
|
200 | << assembleTimer.elapsed() << "(" << 100*assembleTimer.elapsed()/elapsedTot << "%)/" |
| 305 |
2/4✓ Branch 2 taken 100 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 100 times.
✗ Branch 8 not taken.
|
200 | << solveTimer.elapsed() << "(" << 100*solveTimer.elapsed()/elapsedTot << "%)/" |
| 306 |
2/4✓ Branch 2 taken 100 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 100 times.
✗ Branch 8 not taken.
|
200 | << updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)/" |
| 307 |
2/4✓ Branch 2 taken 100 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 100 times.
✗ Branch 8 not taken.
|
200 | << outputTimer.elapsed() << "(" << 100*outputTimer.elapsed()/elapsedTot << "%)" |
| 308 |
1/2✓ Branch 2 taken 100 times.
✗ Branch 3 not taken.
|
200 | << "\n"; |
| 309 | |||
| 310 | // report statistics of this time step | ||
| 311 |
2/4✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
|
200 | timeLoop->reportTimeStep(); |
| 312 | |||
| 313 | // use same time step size in next time step | ||
| 314 | // this function can be used to implement time step adaptivity | ||
| 315 |
2/4✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
|
200 | timeLoop->setTimeStepSize(dt); |
| 316 | } | ||
| 317 | // [[/codeblock]] | ||
| 318 | // | ||
| 319 | // ### Finalize | ||
| 320 | // | ||
| 321 | // [[codeblock]] | ||
| 322 | // end the time loop and produce report | ||
| 323 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
|
2 | timeLoop->finalize(gridView3D.comm()); |
| 324 | |||
| 325 | // output (manual) timings | ||
| 326 |
2/4✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
3 | std::cout << "Transport computation took " << timer.elapsed() - gridSetupTime << " seconds." << std::endl; |
| 327 |
2/4✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
3 | std::cout << "Simulation took " << timer.elapsed() << " seconds." << std::endl; |
| 328 | |||
| 329 | // Print used/unused parameters at the end (good for debugging) | ||
| 330 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | Parameters::print(); |
| 331 | |||
| 332 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
1 | return 0; |
| 333 | } | ||
| 334 | // [[/codeblock]] | ||
| 335 | // [[/content]] | ||
| 336 |