24#ifndef DUMUX_FV_ASSEMBLER_HH
25#define DUMUX_FV_ASSEMBLER_HH
32#include <dune/istl/matrixindexset.hh>
56template<
class DiscretizationMethod>
62 template<
class TypeTag,
class Impl, DiffMethod diffMethod,
bool isImplicit>
69 template<
class TypeTag,
class Impl, DiffMethod diffMethod,
bool isImplicit>
76 template<
class TypeTag,
class Impl, DiffMethod diffMethod,
bool isImplicit>
83 template<
class TypeTag,
class Impl, DiffMethod diffMethod,
bool isImplicit>
90 template<
class TypeTag,
class Impl, DiffMethod diffMethod,
bool isImplicit>
97 template<
class TypeTag,
class Impl, DiffMethod diffMethod,
bool isImplicit>
101template<
class TypeTag,
class Impl, DiffMethod diffMethod,
bool isImplicit>
104>::template type<TypeTag, Impl, diffMethod, isImplicit>;
117template<
class TypeTag, DiffMethod diffMethod,
bool isImplicit = true>
121 using GridView =
typename GridGeo::GridView;
123 using Element =
typename GridView::template Codim<0>::Entity;
124 using ElementSeed =
typename GridView::Grid::template Codim<0>::EntitySeed;
154 , isStationaryProblem_(true)
156 static_assert(isImplicit,
"Explicit assembler for stationary problem doesn't make sense!");
160 && getParam<bool>(
"Assembly.Multithreading",
true);
162 maybeComputeColors_();
173 std::shared_ptr<const TimeLoop> timeLoop,
178 , timeLoop_(timeLoop)
180 , isStationaryProblem_(!timeLoop)
185 && getParam<bool>(
"Assembly.Multithreading",
true);
187 maybeComputeColors_();
194 template<
class PartialReassembler = DefaultPartialReassembler>
197 checkAssemblerState_();
198 resetJacobian_(partialReassembler);
201 assemble_([&](
const Element& element)
203 LocalAssembler localAssembler(*
this, element, curSol);
204 localAssembler.assembleJacobianAndResidual(*jacobian_, *residual_, *gridVariables_, partialReassembler);
207 enforcePeriodicConstraints_(*jacobian_, *residual_, curSol, *gridGeometry_);
215 checkAssemblerState_();
218 assemble_([&](
const Element& element)
220 LocalAssembler localAssembler(*
this, element, curSol);
221 localAssembler.assembleJacobianAndResidual(*jacobian_, *gridVariables_);
235 checkAssemblerState_();
237 assemble_([&](
const Element& element)
239 LocalAssembler localAssembler(*
this, element, curSol);
240 localAssembler.assembleResidual(r);
251 static bool warningIssued =
false;
257 using DM =
typename GridGeometry::VertexMapper;
260 PVHelper vectorHelper(
gridView(), gridGeometry_->vertexMapper());
262 vectorHelper.makeNonOverlappingConsistent(
residual);
265 else if (!warningIssued)
268 std::cout <<
"\nWarning: norm calculation adds entries corresponding to\n"
269 <<
"overlapping entities multiple times. Please use the norm\n"
270 <<
"function provided by a linear solver instead." << std::endl;
272 warningIssued =
true;
278 result2 =
gridView().comm().sum(result2);
281 return sqrt(result2);
290 std::shared_ptr<SolutionVector> r)
296 if (jacobian_->buildMode() == JacobianMatrix::BuildMode::unknown)
297 jacobian_->setBuildMode(JacobianMatrix::random);
298 else if (jacobian_->buildMode() != JacobianMatrix::BuildMode::random)
299 DUNE_THROW(Dune::NotImplemented,
"Only BCRS matrices with random build mode are supported at the moment");
302 setJacobianPattern_();
311 jacobian_ = std::make_shared<JacobianMatrix>();
312 jacobian_->setBuildMode(JacobianMatrix::random);
313 residual_ = std::make_shared<SolutionVector>();
316 setJacobianPattern_();
325 setJacobianPattern_();
326 maybeComputeColors_();
331 {
return gridGeometry_->numDofs(); }
335 {
return *problem_; }
339 {
return *gridGeometry_; }
347 {
return *gridVariables_; }
351 {
return *gridVariables_; }
355 {
return *jacobian_; }
359 {
return *residual_; }
363 {
return *prevSol_; }
370 { timeLoop_ = timeLoop; isStationaryProblem_ = !
static_cast<bool>(timeLoop); }
383 {
return isStationaryProblem_; }
389 {
return LocalResidual(problem_.get(), timeLoop_.get()); }
411 void setJacobianPattern_()
418 const auto occupationPattern = getJacobianPattern<isImplicit>(
gridGeometry());
421 occupationPattern.exportIdx(*jacobian_);
425 void setResidualSize_()
426 { residual_->resize(
numDofs()); }
429 void maybeComputeColors_()
431 if (enableMultithreading_)
436 void resetResidual_()
440 residual_ = std::make_shared<SolutionVector>();
448 template <
class PartialReassembler = DefaultPartialReassembler>
449 void resetJacobian_(
const PartialReassembler *partialReassembler =
nullptr)
453 jacobian_ = std::make_shared<JacobianMatrix>();
454 jacobian_->setBuildMode(JacobianMatrix::random);
455 setJacobianPattern_();
458 if (partialReassembler)
459 partialReassembler->resetJacobian(*
this);
465 void checkAssemblerState_()
const
467 if (!isStationaryProblem_ && !prevSol_)
468 DUNE_THROW(Dune::InvalidStateException,
"Assembling instationary problem but previous solution was not set!");
476 template<
typename AssembleElementFunc>
477 void assemble_(AssembleElementFunc&& assembleElement)
const
480 bool succeeded =
false;
485 if (enableMultithreading_)
487 assert(elementSets_.size() > 0);
493 for (
const auto& elements : elementSets_)
497 const auto element = gridView().grid().entity(elements[i]);
498 assembleElement(element);
503 for (
const auto& element : elements(
gridView()))
504 assembleElement(element);
510 catch (NumericalProblem &e)
512 std::cout <<
"rank " <<
gridView().comm().rank()
513 <<
" caught an exception while assembling:" << e.what()
520 succeeded =
gridView().comm().min(succeeded);
524 DUNE_THROW(NumericalProblem,
"A process did not succeed in linearizing the system");
527 template<
class GG> std::enable_if_t<GG::discMethod == DiscretizationMethods::box, void>
532 if (m.first < m.second)
535 res[m.first] += res[m.second];
536 const auto end = jac[m.second].end();
537 for (
auto it = jac[m.second].begin(); it != end; ++it)
538 jac[m.first][it.index()] += (*it);
541 res[m.second] = curSol[m.second] - curSol[m.first];
542 for (
auto it = jac[m.second].begin(); it != end; ++it)
543 (*it) = it.index() == m.second ? 1.0 : it.index() == m.first ? -1.0 : 0.0;
548 template<
class GG> std::enable_if_t<GG::discMethod != DiscretizationMethods::box, void>
549 enforcePeriodicConstraints_(
JacobianMatrix& jac, SolutionVector& res,
const SolutionVector& curSol,
const GG&
gridGeometry) {}
552 std::shared_ptr<const Problem> problem_;
555 std::shared_ptr<const GridGeometry> gridGeometry_;
558 std::shared_ptr<GridVariables> gridVariables_;
561 std::shared_ptr<const TimeLoop> timeLoop_;
564 const SolutionVector* prevSol_ =
nullptr;
567 bool isStationaryProblem_;
570 std::shared_ptr<JacobianMatrix> jacobian_;
571 std::shared_ptr<SolutionVector> residual_;
574 bool enableMultithreading_ =
false;
575 std::deque<std::vector<ElementSeed>> elementSets_;
An assembler for Jacobian and residual contribution per element (face-centered staggered methods)
An enum class to define various differentiation methods available in order to compute the derivatives...
An assembler for Jacobian and residual contribution per element (face-centered diamond methods)
An assembler for Jacobian and residual contribution per element (cell-centered methods)
Helper function to generate Jacobian pattern for different discretization methods.
An assembler for Jacobian and residual contribution per element (pq1bubble method)
Coloring schemes for shared-memory-parallel assembly.
An assembler for Jacobian and residual contribution per element (box method)
Parallel for loop (multithreading)
dune-grid capabilities compatibility layer
The available discretization methods in Dumux.
Provides a helper class for nonoverlapping decomposition.
constexpr bool isSerial()
Checking whether the backend is serial.
Definition: multithreading.hh:57
void parallelFor(const std::size_t count, const FunctorType &functor)
A parallel for loop (multithreading)
Definition: parallel_for.hh:172
Adaption of the non-isothermal two-phase two-component flow model to problems with CO2.
Definition: adapt.hh:29
typename GetProp< TypeTag, Property >::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:180
auto computeColoring(const GridGeometry &gg, int verbosity=1)
Compute iterable lists of element seeds partitioned by color.
Definition: coloring.hh:251
Distance implementation details.
Definition: cvfelocalresidual.hh:37
typename LocalAssemblerChooser< typename GetPropType< TypeTag, Properties::GridGeometry >::DiscretizationMethod >::template type< TypeTag, Impl, diffMethod, isImplicit > LocalAssemblerChooser_t
Definition: assembly/fvassembler.hh:104
bool supportsMultithreading(const GridView &gridView)
Definition: gridcapabilities.hh:86
CVFE< CVFEMethods::CR_RT > FCDiamond
Definition: method.hh:90
constexpr Box box
Definition: method.hh:136
CVFE< CVFEMethods::PQ1 > Box
Definition: method.hh:83
CVFE< CVFEMethods::PQ1Bubble > PQ1Bubble
Definition: method.hh:97
An assembler for Jacobian and residual contribution per element (box methods)
Definition: boxlocalassembler.hh:257
An assembler for Jacobian and residual contribution per element (cell-centered methods)
Definition: cclocalassembler.hh:136
Traits specifying if a given discretization tag supports coloring.
Definition: coloring.hh:304
An assembler for Jacobian and residual contribution per element (Face-centered methods)
Definition: fcdiamondlocalassembler.hh:290
An assembler for Jacobian and residual contribution per element (Face-centered methods)
Definition: fclocalassembler.hh:296
Definition: assembly/fvassembler.hh:57
A linear system assembler (residual and Jacobian) for finite volume schemes (box, tpfa,...
Definition: assembly/fvassembler.hh:119
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: assembly/fvassembler.hh:134
void assembleJacobian(const SolutionVector &curSol)
Assembles only the global Jacobian of the residual.
Definition: assembly/fvassembler.hh:213
void updateAfterGridAdaption()
Resizes jacobian and residual and recomputes colors.
Definition: assembly/fvassembler.hh:322
const GridView & gridView() const
The gridview.
Definition: assembly/fvassembler.hh:342
const Problem & problem() const
The problem.
Definition: assembly/fvassembler.hh:334
void setLinearSystem()
The version without arguments uses the default constructor to create the jacobian and residual object...
Definition: assembly/fvassembler.hh:309
SolutionVector ResidualType
Definition: assembly/fvassembler.hh:140
bool isStationaryProblem() const
Whether we are assembling a stationary or instationary problem.
Definition: assembly/fvassembler.hh:382
GetPropType< TypeTag, Properties::JacobianMatrix > JacobianMatrix
Definition: assembly/fvassembler.hh:135
LocalResidual localResidual() const
Create a local residual object (used by the local assembler)
Definition: assembly/fvassembler.hh:388
FVAssembler(std::shared_ptr< const Problem > problem, std::shared_ptr< const GridGeometry > gridGeometry, std::shared_ptr< GridVariables > gridVariables, std::shared_ptr< const TimeLoop > timeLoop, const SolutionVector &prevSol)
The constructor for instationary problems.
Definition: assembly/fvassembler.hh:170
Scalar residualNorm(const SolutionVector &curSol) const
compute the residual and return it's vector norm
Definition: assembly/fvassembler.hh:245
std::size_t numDofs() const
Returns the number of degrees of freedom.
Definition: assembly/fvassembler.hh:330
GridVariables & gridVariables()
The global grid variables.
Definition: assembly/fvassembler.hh:346
JacobianMatrix & jacobian()
The jacobian matrix.
Definition: assembly/fvassembler.hh:354
void assembleResidual(ResidualType &r, const SolutionVector &curSol) const
assemble a residual r
Definition: assembly/fvassembler.hh:233
GetPropType< TypeTag, Properties::Problem > Problem
Definition: assembly/fvassembler.hh:137
void setTimeLoop(std::shared_ptr< const TimeLoop > timeLoop)
Set time loop for instationary problems.
Definition: assembly/fvassembler.hh:369
void assembleResidual(const SolutionVector &curSol)
compute the residuals using the internal residual
Definition: assembly/fvassembler.hh:226
GetPropType< TypeTag, Properties::GridVariables > GridVariables
Definition: assembly/fvassembler.hh:138
void setPreviousSolution(const SolutionVector &u)
Sets the solution from which to start the time integration. Has to be called prior to assembly for ti...
Definition: assembly/fvassembler.hh:376
const GridGeometry & gridGeometry() const
The global finite volume geometry.
Definition: assembly/fvassembler.hh:338
void resetTimeStep(const SolutionVector &cursol)
Reset the gridVariables.
Definition: assembly/fvassembler.hh:402
const GridVariables & gridVariables() const
The global grid variables.
Definition: assembly/fvassembler.hh:350
void setLinearSystem(std::shared_ptr< JacobianMatrix > A, std::shared_ptr< SolutionVector > r)
Tells the assembler which jacobian and residual to use. This also resizes the containers to the requi...
Definition: assembly/fvassembler.hh:289
void assembleJacobianAndResidual(const SolutionVector &curSol, const PartialReassembler *partialReassembler=nullptr)
Assembles the global Jacobian of the residual and the residual for the current solution.
Definition: assembly/fvassembler.hh:195
GridGeo GridGeometry
Definition: assembly/fvassembler.hh:136
FVAssembler(std::shared_ptr< const Problem > problem, std::shared_ptr< const GridGeometry > gridGeometry, std::shared_ptr< GridVariables > gridVariables)
The constructor for stationary problems.
Definition: assembly/fvassembler.hh:147
SolutionVector & residual()
The residual vector (rhs)
Definition: assembly/fvassembler.hh:358
const SolutionVector & prevSol() const
The solution of the previous time step.
Definition: assembly/fvassembler.hh:362
void updateGridVariables(const SolutionVector &cursol)
Update the grid variables.
Definition: assembly/fvassembler.hh:394
detects which entries in the Jacobian have to be recomputed
Definition: partialreassembler.hh:432
An assembler for Jacobian and residual contribution per element (PQ1Bubble methods)
Definition: pq1bubblelocalassembler.hh:304
Manages the handling of time dependent problems.
Definition: common/timeloop.hh:68
Definition: parallelhelpers.hh:473
Declares all properties used in Dumux.
Manages the handling of time dependent problems.