94 using Toolbox = MathToolbox<Evaluation>;
96 using Element =
typename GridView::template Codim<0>::Entity;
97 using ElementIterator =
typename GridView::template Codim<0>::Iterator;
99 using Vector = GlobalEqVector;
101 using IstlMatrix =
typename SparseMatrixAdapter::IstlMatrix;
106 using MatrixBlock =
typename SparseMatrixAdapter::MatrixBlock;
107 using VectorBlock = Dune::FieldVector<Scalar, numEq>;
124 auto it = elementCtx_.begin();
125 const auto& endIt = elementCtx_.end();
126 for (; it != endIt; ++it)
147 simulatorPtr_ = &simulator;
149 auto it = elementCtx_.begin();
150 const auto& endIt = elementCtx_.end();
151 for (; it != endIt; ++it){
154 elementCtx_.resize(0);
155 fullDomain_ = std::make_unique<FullDomain>(simulator.
gridView());
200 template <
class SubDomainType>
208 initFirstIteration_();
211 if (
static_cast<std::size_t
>(domain.view.size(0)) == model_().numTotalDof()) {
215 resetSystem_(domain);
223 catch (
const std::exception& e)
225 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
226 <<
" caught an exception while linearizing:" << e.what()
227 <<
"\n" << std::flush;
232 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
233 <<
" caught an exception while linearizing"
234 <<
"\n" << std::flush;
237 succeeded = simulator_().
gridView().comm().min(succeeded);
240 throw NumericalProblem(
"A process did not succeed in linearizing the system");
244 { jacobian_->finalize(); }
256 auto& model = model_();
257 const auto& comm = simulator_().
gridView().comm();
258 for (
unsigned auxModIdx = 0; auxModIdx < model.numAuxiliaryModules(); ++auxModIdx) {
259 bool succeeded =
true;
261 model.auxiliaryModule(auxModIdx)->linearize(*jacobian_, residual_);
263 catch (
const std::exception& e) {
266 std::cout <<
"rank " << simulator_().
gridView().comm().rank()
267 <<
" caught an exception while linearizing:" << e.what()
268 <<
"\n" << std::flush;
271 succeeded = comm.min(succeeded);
274 throw NumericalProblem(
"linearization of an auxiliary equation failed");
282 {
return *jacobian_; }
285 {
return *jacobian_; }
291 {
return residual_; }
294 {
return residual_; }
296 void setLinearizationType(LinearizationType linearizationType){
297 linearizationType_ = linearizationType;
300 const LinearizationType& getLinearizationType()
const{
301 return linearizationType_;
304 void updateDiscretizationParameters()
309 void updateBoundaryConditionData()
314 void updateFlowsInfo() {
324 {
return constraintsMap_; }
340 {
return floresInfo_;}
342 template <
class SubDomainType>
343 void resetSystem_(
const SubDomainType& domain)
346 initFirstIteration_();
350 using GridViewType =
decltype(domain.view);
351 ThreadedEntityIterator<GridViewType, 0> threadedElemIt(domain.view);
357 auto elemIt = threadedElemIt.beginParallel();
358 MatrixBlock zeroBlock;
360 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
361 const Element& elem = *elemIt;
362 ElementContext& elemCtx = *elementCtx_[threadId];
363 elemCtx.updatePrimaryStencil(elem);
365 for (
unsigned primaryDofIdx = 0;
366 primaryDofIdx < elemCtx.numPrimaryDof(0);
369 unsigned globI = elemCtx.globalSpaceIndex(primaryDofIdx, 0);
370 residual_[globI] = 0.0;
371 jacobian_->clearRow(globI, 0.0);
378 Simulator& simulator_()
379 {
return *simulatorPtr_; }
380 const Simulator& simulator_()
const
381 {
return *simulatorPtr_; }
384 {
return simulator_().
problem(); }
385 const Problem& problem_()
const
386 {
return simulator_().
problem(); }
389 {
return simulator_().
model(); }
390 const Model& model_()
const
391 {
return simulator_().
model(); }
393 const GridView& gridView_()
const
394 {
return problem_().gridView(); }
396 const ElementMapper& elementMapper_()
const
397 {
return model_().elementMapper(); }
399 const DofMapper& dofMapper_()
const
400 {
return model_().dofMapper(); }
402 void initFirstIteration_()
408 residual_.resize(model_().numTotalDof());
414 elementCtx_[threadId] =
new ElementContext(simulator_());
420 const auto& model = model_();
421 Stencil stencil(gridView_(), model_().dofMapper());
425 sparsityPattern_.clear();
426 sparsityPattern_.resize(model.numTotalDof());
428 for (
const auto& elem : elements(gridView_())) {
429 stencil.update(elem);
431 for (
unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
432 unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
434 for (
unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
435 unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
436 sparsityPattern_[myIdx].insert(neighborIdx);
443 size_t numAuxMod = model.numAuxiliaryModules();
444 for (
unsigned auxModIdx = 0; auxModIdx < numAuxMod; ++auxModIdx)
445 model.auxiliaryModule(auxModIdx)->addNeighbors(sparsityPattern_);
448 jacobian_.reset(
new SparseMatrixAdapter(simulator_()));
451 jacobian_->reserve(sparsityPattern_);
464 void updateConstraintsMap_()
466 if (!enableConstraints_())
470 constraintsMap_.clear();
473 ThreadedEntityIterator<GridView, 0> threadedElemIt(gridView_());
479 ElementIterator elemIt = threadedElemIt.beginParallel();
480 for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
483 const Element& elem = *elemIt;
484 ElementContext& elemCtx = *elementCtx_[threadId];
485 elemCtx.updateStencil(elem);
489 for (
unsigned primaryDofIdx = 0;
490 primaryDofIdx < elemCtx.numPrimaryDof(0);
493 Constraints constraints;
494 elemCtx.problem().constraints(constraints,
498 if (constraints.isActive()) {
499 unsigned globI = elemCtx.globalSpaceIndex(primaryDofIdx, 0);
500 constraintsMap_[globI] = constraints;
509 template <
class SubDomainType>
510 void linearize_(
const SubDomainType& domain)
512 OPM_TIMEBLOCK(linearize_);
521 if (model_().newtonMethod().numIterations() == 0)
522 updateConstraintsMap_();
524 applyConstraintsToSolution_();
529 std::mutex exceptionLock;
533 std::exception_ptr exceptionPtr =
nullptr;
536 using GridViewType =
decltype(domain.view);
537 ThreadedEntityIterator<GridViewType, 0> threadedElemIt(domain.view);
542 auto elemIt = threadedElemIt.beginParallel();
543 auto nextElemIt = elemIt;
545 for (; !threadedElemIt.isFinished(elemIt); elemIt = nextElemIt) {
548 nextElemIt = threadedElemIt.increment();
549 if (!threadedElemIt.isFinished(nextElemIt)) {
550 const auto& nextElem = *nextElemIt;
551 if (linearizeNonLocalElements
552 || nextElem.partitionType() == Dune::InteriorEntity)
554 model_().prefetch(nextElem);
555 problem_().prefetch(nextElem);
559 const auto& elem = *elemIt;
560 if (!linearizeNonLocalElements && elem.partitionType() != Dune::InteriorEntity)
563 linearizeElement_(elem);
576 std::lock_guard<std::mutex> take(exceptionLock);
577 exceptionPtr = std::current_exception();
578 threadedElemIt.setFinished();
586 std::rethrow_exception(exceptionPtr);
589 applyConstraintsToLinearization_();
594 template <
class ElementType>
599 ElementContext *elementCtx = elementCtx_[threadId];
600 auto& localLinearizer = model_().localLinearizer(threadId);
603 localLinearizer.linearize(*elementCtx, elem);
607 globalMatrixMutex_.lock();
609 size_t numPrimaryDof = elementCtx->numPrimaryDof(0);
610 for (
unsigned primaryDofIdx = 0; primaryDofIdx < numPrimaryDof; ++ primaryDofIdx) {
611 unsigned globI = elementCtx->globalSpaceIndex(primaryDofIdx, 0);
614 residual_[globI] += localLinearizer.residual(primaryDofIdx);
617 for (
unsigned dofIdx = 0; dofIdx < elementCtx->numDof(0); ++ dofIdx) {
618 unsigned globJ = elementCtx->globalSpaceIndex(dofIdx, 0);
620 jacobian_->addToBlock(globJ, globI, localLinearizer.jacobian(dofIdx, primaryDofIdx));
625 globalMatrixMutex_.unlock();
630 void applyConstraintsToSolution_()
632 if (!enableConstraints_())
636 auto& sol = model_().solution(0);
637 auto& oldSol = model_().solution(1);
639 auto it = constraintsMap_.begin();
640 const auto& endIt = constraintsMap_.end();
641 for (; it != endIt; ++it) {
642 sol[it->first] = it->second;
643 oldSol[it->first] = it->second;
649 void applyConstraintsToLinearization_()
651 if (!enableConstraints_())
654 auto it = constraintsMap_.begin();
655 const auto& endIt = constraintsMap_.end();
656 for (; it != endIt; ++it) {
657 unsigned constraintDofIdx = it->first;
661 jacobian_->clearRow(constraintDofIdx, Scalar(1.0));
664 residual_[constraintDofIdx] = 0.0;
668 static bool enableConstraints_()
671 Simulator *simulatorPtr_;
672 std::vector<ElementContext*> elementCtx_;
676 std::map<unsigned, Constraints> constraintsMap_;
685 SparseTable<FlowInfo> flowsInfo_;
686 SparseTable<FlowInfo> floresInfo_;
689 std::unique_ptr<SparseMatrixAdapter> jacobian_;
692 GlobalEqVector residual_;
694 LinearizationType linearizationType_;
696 std::mutex globalMatrixMutex_;
698 std::vector<std::set<unsigned int>> sparsityPattern_;
702 explicit FullDomain(
const GridView& v) : view (v) {}
704 std::vector<bool> interior;
709 std::unique_ptr<FullDomain> fullDomain_;