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()) {
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;
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();
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;
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>
346 initFirstIteration_();
362 ElementContext& elemCtx = *elementCtx_[threadId];
363 elemCtx.updatePrimaryStencil(
elem);
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());
429 stencil.update(
elem);
434 for (
unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
435 unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
443 size_t numAuxMod = model.numAuxiliaryModules();
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();
484 ElementContext& elemCtx = *elementCtx_[threadId];
485 elemCtx.updateStencil(
elem);
493 Constraints constraints;
494 elemCtx.problem().constraints(constraints,
498 if (constraints.isActive()) {
500 constraintsMap_[
globI] = constraints;
509 template <
class SubDomainType>
521 if (model_().newtonMethod().numIterations() == 0)
522 updateConstraintsMap_();
524 applyConstraintsToSolution_();
552 ||
nextElem.partitionType() == Dune::InteriorEntity)
563 linearizeElement_(
elem);
589 applyConstraintsToLinearization_();
594 template <
class ElementType>
599 ElementContext *
elementCtx = elementCtx_[threadId];
600 auto& localLinearizer = model_().localLinearizer(threadId);
607 globalMatrixMutex_.lock();
609 size_t numPrimaryDof =
elementCtx->numPrimaryDof(0);
617 for (
unsigned dofIdx = 0; dofIdx <
elementCtx->numDof(0); ++ dofIdx) {
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) {
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_;