From d39a38f86a25d9685fcb2114924ae6f33e1d44cd Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 29 Aug 2019 14:11:20 +0200 Subject: [PATCH 01/10] [multidomain][cc] Rename ImplInverse to Impl. There is only one implementation. --- dumux/multidomain/subdomaincclocalassembler.hh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/dumux/multidomain/subdomaincclocalassembler.hh b/dumux/multidomain/subdomaincclocalassembler.hh index 12bd85cde5..f9b3ce0897 100644 --- a/dumux/multidomain/subdomaincclocalassembler.hh +++ b/dumux/multidomain/subdomaincclocalassembler.hh @@ -120,7 +120,7 @@ public: // for the diagonal jacobian block const auto globalI = this->fvGeometry().gridGeometry().elementMapper().index(this->element()); - res[globalI] = this->asImp_().assembleJacobianAndResidualImplInverse(jacRow[domainId], *std::get(gridVariables)); + res[globalI] = this->asImp_().assembleJacobianAndResidualImpl(jacRow[domainId], *std::get(gridVariables)); // for the coupling blocks using namespace Dune::Hybrid; @@ -308,7 +308,7 @@ public: * \return The element residual at the current solution. */ template - LocalResidualValues assembleJacobianAndResidualImplInverse(JacobianMatrixDiagBlock& A, GridVariables& gridVariables) + LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrixDiagBlock& A, GridVariables& gridVariables) { ////////////////////////////////////////////////////////////////////////////////////////////////// // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the // @@ -640,7 +640,7 @@ public: * \return The element residual at the current solution. */ template - LocalResidualValues assembleJacobianAndResidualImplInverse(JacobianMatrixDiagBlock& A, GridVariables& gridVariables) + LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrixDiagBlock& A, GridVariables& gridVariables) { if (this->assembler().isStationaryProblem()) DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual"); -- GitLab From f2a2a5e25a3c2a97d6673b1be43c164c98c3e9d9 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 29 Aug 2019 15:30:12 +0200 Subject: [PATCH 02/10] [multidomain][couplingmanager] Add throwing default impl of addCouplingDerivatives --- dumux/multidomain/couplingmanager.hh | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/dumux/multidomain/couplingmanager.hh b/dumux/multidomain/couplingmanager.hh index b1fbbb6296..d6874c0f8d 100644 --- a/dumux/multidomain/couplingmanager.hh +++ b/dumux/multidomain/couplingmanager.hh @@ -239,6 +239,28 @@ public: return NumericEpsilon(paramGroup); } + /*! + * \ingroup MultiDomain + * \brief evaluates the element residual derivative of a coupled element of domain i + * with respect to all coupled degrees of freedom of domain j + * + * \param Aij the coupling block of Jacobian between domain i and j + * \param domainI the domain index of domain i + * \param localAssemblerI the local assembler assembling the element residual of an element of domain i + * \param domainJ the domain index of domain j + * + * \note The element whose residual derivatives are to be evaluated can be retrieved from the local assembler + * as localAssemblerI.element() as well as all up-to-date variables and caches. + */ + template + void addCouplingDerivatives(JacobianCouplingBlock& Aij, + Dune::index_constant domainI, + const LocalAssemblerI& localAssemblerI, + Dune::index_constant domainJ) const + { + DUNE_THROW(Dune::NotImplemented, "This coupling manager does not implement analytic derivatives. Use numeric differentiation instead."); + } + /*! * \brief set the pointers to the sub problems * \param problems A tuple of shared pointers to the sub problems -- GitLab From f3e1605bbd490b3739cb63a970fc4dadacee9000 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 29 Aug 2019 15:30:58 +0200 Subject: [PATCH 03/10] [multidomain][box][cc] Implement assembly routine for analytic Jacobians --- .../multidomain/subdomainboxlocalassembler.hh | 105 ++++++++++++++++++ .../multidomain/subdomaincclocalassembler.hh | 36 ++---- 2 files changed, 117 insertions(+), 24 deletions(-) diff --git a/dumux/multidomain/subdomainboxlocalassembler.hh b/dumux/multidomain/subdomainboxlocalassembler.hh index e70259f712..6b7fea2f31 100644 --- a/dumux/multidomain/subdomainboxlocalassembler.hh +++ b/dumux/multidomain/subdomainboxlocalassembler.hh @@ -693,6 +693,111 @@ public: {} }; +/*! + * \ingroup Assembly + * \ingroup BoxDiscretization + * \brief Box local assembler using analytic differentiation and implicit time discretization + */ +template +class SubDomainBoxLocalAssembler +: public SubDomainBoxLocalAssemblerBase, true> +{ + using ThisType = SubDomainBoxLocalAssembler; + using ParentType = SubDomainBoxLocalAssemblerBase; + using GridVariables = GetPropType; + using JacobianMatrix = GetPropType; + using LocalResidual = GetPropType; + using ElementResidualVector = typename LocalResidual::ElementResidualVector; + static constexpr auto domainI = Dune::index_constant(); + +public: + using ParentType::ParentType; + + /*! + * \brief Computes the derivatives with respect to the given element and adds them + * to the global matrix. + * + * \return The element residual at the current solution. + */ + ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables) + { + // get some aliases for convenience + const auto& element = this->element(); + const auto& fvGeometry = this->fvGeometry(); + const auto& problem = this->problem(); + const auto& curElemVolVars = this->curElemVolVars(); + const auto& elemFluxVarsCache = this->elemFluxVarsCache(); + + // get the vecor of the acutal element residuals + const auto origResiduals = this->evalLocalResidual(); + + ////////////////////////////////////////////////////////////////////////////////////////////////// + // // + // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the // + // neighboring elements we do so by computing the derivatives of the fluxes which depend on the // + // actual element. In the actual element we evaluate the derivative of the entire residual. // + // // + ////////////////////////////////////////////////////////////////////////////////////////////////// + + // calculation of the source and storage derivatives + for (const auto& scv : scvs(fvGeometry)) + { + // dof index and corresponding actual pri vars + const auto dofIdx = scv.dofIndex(); + const auto& volVars = curElemVolVars[scv]; + + // derivative of this scv residual w.r.t the d.o.f. of the same scv (because of mass lumping) + // only if the problem is instationary we add derivative of storage term + // TODO if e.g. porosity depends on all dofs in the element, we would have off-diagonal matrix entries!? + if (!this->assembler().isStationaryProblem()) + this->localResidual().addStorageDerivatives(A[dofIdx][dofIdx], problem, element, fvGeometry, volVars, scv); + + // derivative of this scv residual w.r.t the d.o.f. of the same scv (because of mass lumping) + // add source term derivatives + this->localResidual().addSourceDerivatives(A[dofIdx][dofIdx], problem, element, fvGeometry, volVars, scv); + } + + // localJacobian[scvIdx][otherScvIdx][eqIdx][priVarIdx] of the fluxes + for (const auto& scvf : scvfs(fvGeometry)) + { + // add flux term derivatives + if (!scvf.boundary()) + this->localResidual().addFluxDerivatives(A, problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf); + + // the boundary gets special treatment to simplify + // for the user + else + { + // add flux term derivatives + const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); + if (this->elemBcTypes()[insideScv.localDofIndex()].hasNeumann()) + this->localResidual().addRobinFluxDerivatives(A[insideScv.dofIndex()], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf); + } + } + + return origResiduals; + } + + /*! + * \brief Computes the derivatives with respect to the given element and adds them + * to the global matrix. + * + * \return The element residual at the current solution. + */ + template + void assembleJacobianCoupling(Dune::index_constant domainJ, JacobianBlock& A, + const ElementResidualVector& res, GridVariables& gridVariables) + { + /////////////////////////////////////////////////////// + // Calculate derivatives of all dofs in the element // + // with respect to all dofs in the coupling stencil. // + /////////////////////////////////////////////////////// + this->couplingManager().addCouplingDerivatives(A, domainI, *this, domainJ); + } + +}; // implicit assebler with analytic Jacobian + } // end namespace Dumux #endif diff --git a/dumux/multidomain/subdomaincclocalassembler.hh b/dumux/multidomain/subdomaincclocalassembler.hh index f9b3ce0897..79783b5b1a 100644 --- a/dumux/multidomain/subdomaincclocalassembler.hh +++ b/dumux/multidomain/subdomaincclocalassembler.hh @@ -771,8 +771,8 @@ class SubDomainCCLocalAssembler::Entity; using Problem = GetPropType; - enum { numEq = GetPropType::numEq() }; - enum { dim = GridView::dimension }; + static constexpr int numEq = GetPropType::numEq(); + static constexpr int dim = GridView::dimension; static constexpr auto domainI = Dune::index_constant(); @@ -873,42 +873,30 @@ public: void assembleJacobianCoupling(Dune::index_constant domainJ, JacobianBlock& A, const LocalResidualValues& res, GridVariables& gridVariables) { - //////////////////////////////////////////////////////////////////////////////////////////////////////// - // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. // - //////////////////////////////////////////////////////////////////////////////////////////////////////// - - // get some aliases for convenience - const auto& element = this->element(); - const auto& fvGeometry = this->fvGeometry(); - const auto& gridGeometry = fvGeometry.gridGeometry(); - auto&& curElemVolVars = this->curElemVolVars(); - // auto&& elemFluxVarsCache = this->elemFluxVarsCache(); + /////////////////////////////////////////////////////// + // Calculate derivatives of all dofs in the element // + // with respect to all dofs in the coupling stencil. // + /////////////////////////////////////////////////////// + this->couplingManager().addCouplingDerivatives(A, domainI, *this, domainJ); - // get stencil informations - const auto globalI = gridGeometry.elementMapper().index(element); - const auto& stencil = this->couplingManager().couplingStencil(domainI, element, domainJ); - - for (const auto globalJ : stencil) + // add the current partial derivatives to the global jacobian matrix + if constexpr (Problem::enableInternalDirichletConstraints()) { - const auto& elementJ = this->assembler().gridGeometry(domainJ).element(globalJ); - this->couplingManager().addCouplingDerivatives(A[globalI][globalJ], domainI, element, fvGeometry, curElemVolVars, domainJ, elementJ); - // add the current partial derivatives to the global jacobian matrix - if constexpr (Problem::enableInternalDirichletConstraints()) + for (const auto globalJ : stencil) { const auto& scv = fvGeometry.scv(globalI); const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv); if (internalDirichletConstraints.any()) - { for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) if (internalDirichletConstraints[eqIdx]) A[globalI][globalJ][eqIdx][pvIdx] = 0.0; - } } } } -}; + +}; // implicit assebler with analytic Jacobian } // end namespace Dumux -- GitLab From 7420016e509a2f961e2ccacea1d51b3af5ae48d3 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 29 Aug 2019 20:36:01 +0200 Subject: [PATCH 04/10] [couplingmanager][doc] Improve docstring --- dumux/multidomain/couplingmanager.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dumux/multidomain/couplingmanager.hh b/dumux/multidomain/couplingmanager.hh index d6874c0f8d..4583039bc0 100644 --- a/dumux/multidomain/couplingmanager.hh +++ b/dumux/multidomain/couplingmanager.hh @@ -241,8 +241,8 @@ public: /*! * \ingroup MultiDomain - * \brief evaluates the element residual derivative of a coupled element of domain i - * with respect to all coupled degrees of freedom of domain j + * \brief Evaluates the element residual derivative of a coupled element of domain i + * with respect to all coupled degrees of freedom of domain j (for analytic Jacobians) * * \param Aij the coupling block of Jacobian between domain i and j * \param domainI the domain index of domain i -- GitLab From daf1a605e27693e6f7af0e15a1e52dbefbca83df Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 29 Aug 2019 20:36:50 +0200 Subject: [PATCH 05/10] [assembly] Pass element volume variables to source derivative function --- dumux/assembly/boxlocalassembler.hh | 2 +- dumux/assembly/cclocalassembler.hh | 2 +- dumux/common/fvproblem.hh | 2 +- dumux/porousmediumflow/1p/incompressiblelocalresidual.hh | 4 ++-- dumux/porousmediumflow/2p/incompressiblelocalresidual.hh | 2 +- dumux/porousmediumflow/richards/localresidual.hh | 6 ++++-- dumux/porousmediumflow/tracer/localresidual.hh | 2 +- 7 files changed, 11 insertions(+), 9 deletions(-) diff --git a/dumux/assembly/boxlocalassembler.hh b/dumux/assembly/boxlocalassembler.hh index 292ddf3be6..7cec280faf 100644 --- a/dumux/assembly/boxlocalassembler.hh +++ b/dumux/assembly/boxlocalassembler.hh @@ -560,7 +560,7 @@ public: problem, element, fvGeometry, - volVars, + curElemVolVars, scv); } diff --git a/dumux/assembly/cclocalassembler.hh b/dumux/assembly/cclocalassembler.hh index 9a9b8168e2..6855b933e8 100644 --- a/dumux/assembly/cclocalassembler.hh +++ b/dumux/assembly/cclocalassembler.hh @@ -540,7 +540,7 @@ public: this->localResidual().addStorageDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv); // add source term derivatives - this->localResidual().addSourceDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv); + this->localResidual().addSourceDerivatives(A[globalI][globalI], problem, element, fvGeometry, curElemVolVars, scv); // add flux derivatives for each scvf for (const auto& scvf : scvfs(fvGeometry)) diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh index f306356947..6b27cc5555 100644 --- a/dumux/common/fvproblem.hh +++ b/dumux/common/fvproblem.hh @@ -419,7 +419,7 @@ public: void addSourceDerivatives(MatrixBlock& block, const Element& element, const FVElementGeometry& fvGeometry, - const VolumeVariables& volVars, + const ElementVolumeVariables& curElemVolVars, const SubControlVolume& scv) const {} /*! diff --git a/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh b/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh index 184fffe068..7e016d799c 100644 --- a/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh +++ b/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh @@ -76,10 +76,10 @@ public: const Problem& problem, const Element& element, const FVElementGeometry& fvGeometry, - const VolumeVariables& curVolVars, + const ElementVolumeVariables& curElemVolVars, const SubControlVolume& scv) const { - problem.addSourceDerivatives(partialDerivatives, element, fvGeometry, curVolVars, scv); + problem.addSourceDerivatives(partialDerivatives, element, fvGeometry, curElemVolVars, scv); } //! Flux derivatives for the cell-centered tpfa scheme diff --git a/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh b/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh index 9475f5c797..f8658eb5cd 100644 --- a/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh +++ b/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh @@ -128,7 +128,7 @@ public: const Problem& problem, const Element& element, const FVElementGeometry& fvGeometry, - const VolumeVariables& curVolVars, + const ElementVolumeVariables& curElemVolVars, const SubControlVolume& scv) const { /* TODO maybe forward to problem for the user to implement the source derivatives?*/ } diff --git a/dumux/porousmediumflow/richards/localresidual.hh b/dumux/porousmediumflow/richards/localresidual.hh index 5434833c98..31a121ae2d 100644 --- a/dumux/porousmediumflow/richards/localresidual.hh +++ b/dumux/porousmediumflow/richards/localresidual.hh @@ -224,9 +224,11 @@ public: const Problem& problem, const Element& element, const FVElementGeometry& fvGeometry, - const VolumeVariables& curVolVars, + const ElementVolumeVariables& curElemVolVars, const SubControlVolume& scv) const - { /* TODO maybe forward to problem for the user to implement the source derivatives?*/ } + { + problem.addSourceDerivatives(partialDerivatives, element, fvGeometry, curElemVolVars, scv); + } /*! * \brief Adds flux derivatives for wetting and non-wetting phase for cell-centered FVM using TPFA diff --git a/dumux/porousmediumflow/tracer/localresidual.hh b/dumux/porousmediumflow/tracer/localresidual.hh index 702ff39e99..047ff12eff 100644 --- a/dumux/porousmediumflow/tracer/localresidual.hh +++ b/dumux/porousmediumflow/tracer/localresidual.hh @@ -229,7 +229,7 @@ public: const Problem& problem, const Element& element, const FVElementGeometry& fvGeometry, - const VolumeVariables& curVolVars, + const ElementVolumeVariables& curElemVolVars, const SubControlVolume& scv) const { // TODO maybe forward to the problem? -> necessary for reaction terms -- GitLab From 35aebe8d66e9f415dfe03e082088f7fb5b1e6fdb Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 29 Aug 2019 20:37:32 +0200 Subject: [PATCH 06/10] [multidomain] Enable replacing analytic coupling and source derivative by numeric differences --- .../multidomain/subdomainboxlocalassembler.hh | 195 ++++++++++++++++- .../multidomain/subdomaincclocalassembler.hh | 198 +++++++++++++++++- 2 files changed, 382 insertions(+), 11 deletions(-) diff --git a/dumux/multidomain/subdomainboxlocalassembler.hh b/dumux/multidomain/subdomainboxlocalassembler.hh index 6b7fea2f31..df9e7bbbac 100644 --- a/dumux/multidomain/subdomainboxlocalassembler.hh +++ b/dumux/multidomain/subdomainboxlocalassembler.hh @@ -705,11 +705,15 @@ class SubDomainBoxLocalAssembler; using ParentType = SubDomainBoxLocalAssemblerBase; + using Scalar = GetPropType; using GridVariables = GetPropType; using JacobianMatrix = GetPropType; - using LocalResidual = GetPropType; - using ElementResidualVector = typename LocalResidual::ElementResidualVector; + using ElementResidualVector = typename GetPropType::ElementResidualVector; static constexpr auto domainI = Dune::index_constant(); + static constexpr int numEq = GetPropType::numEq(); + static constexpr bool enableGridFluxVarsCache = getPropValue(); + static constexpr bool enableGridVolVarsCache = getPropValue(); + enum { dim = GetPropType::dimension }; public: using ParentType::ParentType; @@ -725,12 +729,17 @@ public: // get some aliases for convenience const auto& element = this->element(); const auto& fvGeometry = this->fvGeometry(); + const auto& curSol = this->curSol()[domainI]; const auto& problem = this->problem(); - const auto& curElemVolVars = this->curElemVolVars(); + auto&& curElemVolVars = this->curElemVolVars(); const auto& elemFluxVarsCache = this->elemFluxVarsCache(); // get the vecor of the acutal element residuals const auto origResiduals = this->evalLocalResidual(); + auto origSourceResidual = origResiduals; origSourceResidual = 0.0; + static const bool useNumDiff = getParam("Assembly.UseNumericDifferentiationForSourceDerivatives", false); + if (useNumDiff) + origSourceResidual = this->evalLocalSourceResidual(element, curElemVolVars); ////////////////////////////////////////////////////////////////////////////////////////////////// // // @@ -739,6 +748,12 @@ public: // actual element. In the actual element we evaluate the derivative of the entire residual. // // // ////////////////////////////////////////////////////////////////////////////////////////////////// + // create the element solution + auto elemSol = elementSolution(element, curSol, fvGeometry.fvGridGeometry()); + + // create the vector storing the partial derivatives + auto partialDerivs = origResiduals; + partialDerivs = 0.0; // calculation of the source and storage derivatives for (const auto& scv : scvs(fvGeometry)) @@ -755,7 +770,56 @@ public: // derivative of this scv residual w.r.t the d.o.f. of the same scv (because of mass lumping) // add source term derivatives - this->localResidual().addSourceDerivatives(A[dofIdx][dofIdx], problem, element, fvGeometry, volVars, scv); + if (useNumDiff) + { + auto& curVolVars = this->getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv); + const auto origVolVars = curVolVars; + + // calculate derivatives w.r.t to the privars at the dof at hand + for (int pvIdx = 0; pvIdx < numEq; pvIdx++) + { + partialDerivs = 0.0; + + auto evalSource = [&](Scalar priVar) + { + const auto localDofIndex = scv.localDofIndex(); + elemSol[localDofIndex][pvIdx] = priVar; + curVolVars.update(elemSol, this->problem(), element, scv); + this->couplingManager().updateCouplingContext(domainI, *this, domainI, scv.dofIndex(), elemSol[localDofIndex], pvIdx); + return this->evalLocalSourceResidual(element, curElemVolVars); + }; + + // derive the residuals numerically + static const NumericEpsilon eps_{this->problem().paramGroup()}; + static const int numDiffMethod = getParamFromGroup(this->problem().paramGroup(), "Assembly.NumericDifferenceMethod"); + NumericDifferentiation::partialDerivative(evalSource, elemSol[scv.localDofIndex()][pvIdx], partialDerivs, origSourceResidual, + eps_(elemSol[scv.localDofIndex()][pvIdx], pvIdx), numDiffMethod); + + // update the global stiffness matrix with the current partial derivatives + for (auto&& scvJ : scvs(fvGeometry)) + { + for (int eqIdx = 0; eqIdx < numEq; eqIdx++) + { + // A[i][col][eqIdx][pvIdx] is the rate of change of + // the residual of equation 'eqIdx' at dof 'i' + // depending on the primary variable 'pvIdx' at dof + // 'col'. + A[scvJ.dofIndex()][dofIdx][eqIdx][pvIdx] += partialDerivs[scvJ.localDofIndex()][eqIdx]; + } + } + + // restore the original state of the scv's volume variables + curVolVars = origVolVars; + + // restore the original element solution + elemSol[scv.localDofIndex()][pvIdx] = curSol[scv.dofIndex()][pvIdx]; + + // restore coupling manager + this->couplingManager().updateCouplingContext(domainI, *this, domainI, scv.dofIndex(), elemSol[scv.localDofIndex()], pvIdx); + } + } + else + this->localResidual().addSourceDerivatives(A[dofIdx][dofIdx], problem, element, fvGeometry, curElemVolVars, scv); } // localJacobian[scvIdx][otherScvIdx][eqIdx][priVarIdx] of the fluxes @@ -776,6 +840,9 @@ public: } } + // evaluate additional derivatives that might arise from the coupling + this->couplingManager().evalAdditionalDomainDerivatives(domainI, *this, origResiduals, A, gridVariables); + return origResiduals; } @@ -793,7 +860,125 @@ public: // Calculate derivatives of all dofs in the element // // with respect to all dofs in the coupling stencil. // /////////////////////////////////////////////////////// - this->couplingManager().addCouplingDerivatives(A, domainI, *this, domainJ); + static const bool useNumDiff = getParam("Assembly.UseNumericDifferentiationForCouplingDerivatives", false); + if (useNumDiff) + assembleJacobianCouplingNumeric(domainJ, A, res, gridVariables); + else + this->couplingManager().addCouplingDerivatives(A, domainI, *this, domainJ); + } + + /*! + * \brief Computes the derivatives with respect to the given element and adds them + * to the global matrix. + * \todo This is just a copy of the assembleJacobianCoupling of the numeric-diff local assembler + * \return The element residual at the current solution. + */ + template + void assembleJacobianCouplingNumeric(Dune::index_constant domainJ, JacobianBlock& A, + const ElementResidualVector& res, GridVariables& gridVariables) + { + //////////////////////////////////////////////////////////////////////////////////////////////////////// + // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. // + //////////////////////////////////////////////////////////////////////////////////////////////////////// + + // get some aliases for convenience + const auto& element = this->element(); + const auto& fvGeometry = this->fvGeometry(); + auto&& curElemVolVars = this->curElemVolVars(); + auto&& elemFluxVarsCache = this->elemFluxVarsCache(); + + // get element stencil informations + const auto& stencil = this->couplingManager().couplingStencil(domainI, element, domainJ); + + // convenience lambda for call to update self + auto updateSelf = [&] () + { + // Update ourself after the context has been modified. Depending on the + // type of caching, other objects might have to be updated. All ifs can be optimized away. + if (enableGridFluxVarsCache) + { + if (enableGridVolVarsCache) + this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), gridVariables.gridFluxVarsCache()); + else + this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, gridVariables.gridFluxVarsCache()); + } + else + { + if (enableGridVolVarsCache) + this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), elemFluxVarsCache); + else + this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, elemFluxVarsCache); + } + }; + + const auto& curSolJ = this->curSol()[domainJ]; + for (const auto globalJ : stencil) + { + // undeflected privars and privars to be deflected + const auto origPriVarsJ = curSolJ[globalJ]; + auto priVarsJ = origPriVarsJ; + + // the undeflected coupling residual + const auto origResidual = this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ); + + for (int pvIdx = 0; pvIdx < JacobianBlock::block_type::cols; ++pvIdx) + { + auto evalCouplingResidual = [&](Scalar priVar) + { + priVarsJ[pvIdx] = priVar; + this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx); + updateSelf(); + return this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ); + }; + + // derive the residuals numerically + ElementResidualVector partialDerivs(element.subEntities(dim)); + + const auto& paramGroup = this->assembler().problem(domainJ).paramGroup(); + static const int numDiffMethod = getParamFromGroup(paramGroup, "Assembly.NumericDifferenceMethod"); + static const auto epsCoupl = this->couplingManager().numericEpsilon(domainJ, paramGroup); + + NumericDifferentiation::partialDerivative(evalCouplingResidual, origPriVarsJ[pvIdx], partialDerivs, origResidual, + epsCoupl(origPriVarsJ[pvIdx], pvIdx), numDiffMethod); + + // update the global stiffness matrix with the current partial derivatives + for (auto&& scv : scvs(fvGeometry)) + { + for (int eqIdx = 0; eqIdx < numEq; eqIdx++) + { + // A[i][col][eqIdx][pvIdx] is the rate of change of + // the residual of equation 'eqIdx' at dof 'i' + // depending on the primary variable 'pvIdx' at dof + // 'col'. + + // If the dof is coupled by a Dirichlet condition, + // set the derived value only once (i.e. overwrite existing values). + // For other dofs, add the contribution of the partial derivative. + const auto bcTypes = this->elemBcTypes()[scv.localDofIndex()]; + if (bcTypes.isCouplingDirichlet(eqIdx)) + A[scv.dofIndex()][globalJ][eqIdx][pvIdx] = partialDerivs[scv.localDofIndex()][eqIdx]; + else if (bcTypes.isDirichlet(eqIdx)) + A[scv.dofIndex()][globalJ][eqIdx][pvIdx] = 0.0; + else + A[scv.dofIndex()][globalJ][eqIdx][pvIdx] += partialDerivs[scv.localDofIndex()][eqIdx]; + } + } + + // restore the current element solution + priVarsJ[pvIdx] = origPriVarsJ[pvIdx]; + + // restore the undeflected state of the coupling context + this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx); + } + } + + // restore original state of the flux vars cache and/or vol vars in case of global caching. + // This has to be done in order to guarantee that everything is in an undeflected + // state before the assembly of another element is called. In the case of local caching + // this is obsolete because the local views used here go out of scope after this. + // We only have to do this for the last primary variable, for all others the flux var cache + // is updated with the correct element volume variables before residual evaluations + updateSelf(); } }; // implicit assebler with analytic Jacobian diff --git a/dumux/multidomain/subdomaincclocalassembler.hh b/dumux/multidomain/subdomaincclocalassembler.hh index 79783b5b1a..501dc85639 100644 --- a/dumux/multidomain/subdomaincclocalassembler.hh +++ b/dumux/multidomain/subdomaincclocalassembler.hh @@ -775,6 +775,9 @@ class SubDomainCCLocalAssembler(); + static constexpr int numEq = GetPropType::numEq(); + static constexpr bool enableGridFluxVarsCache = getPropValue(); + static constexpr bool enableGridVolVarsCache = getPropValue(); public: using ParentType::ParentType; @@ -792,11 +795,12 @@ public: const auto& problem = this->problem(); const auto& element = this->element(); const auto& fvGeometry = this->fvGeometry(); - const auto& curElemVolVars = this->curElemVolVars(); + const auto& fvGridGeometry = fvGeometry.fvGridGeometry(); + auto&& curElemVolVars = this->curElemVolVars(); const auto& elemFluxVarsCache = this->elemFluxVarsCache(); // get reference to the element's current vol vars - const auto globalI = this->assembler().gridGeometry(domainI).elementMapper().index(element); + const auto globalI = fvGridGeometry.elementMapper().index(element); const auto& scv = fvGeometry.scv(globalI); const auto& volVars = curElemVolVars[scv]; @@ -805,7 +809,80 @@ public: this->localResidual().addStorageDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv); // add source term derivatives - this->localResidual().addSourceDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv); + static const bool useNumDiff = getParam("Assembly.UseNumericDifferentiationForSourceDerivatives", false); + if (useNumDiff) + { + const auto sourceResidual = this->evalLocalSourceResidual(element, curElemVolVars)[0]; + auto& curVolVars = ParentType::getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv); + + // save a copy of the original privars and vol vars in order + // to restore the original solution after deflection + const auto& curSol = this->curSol()[domainI]; + const auto origPriVars = curSol[globalI]; + const auto origVolVars = curVolVars; + + // element solution container to be deflected + auto elemSol = elementSolution(element, curSol, fvGridGeometry); + + // derivatives in the neighbors with repect to the current elements + LocalResidualValues partialDeriv; + for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) + { + // reset derivatives of element dof with respect to itself + partialDeriv = 0.0; + auto evalSource = [&](Scalar priVar) + { + // update the volume variables and calculate + // the residual with the deflected primary variables + elemSol[0][pvIdx] = priVar; + this->couplingManager().updateCouplingContext(domainI, *this, domainI, globalI, elemSol[0], pvIdx); + curVolVars.update(elemSol, this->problem(), element, scv); + return this->evalLocalSourceResidual(element, curElemVolVars)[0]; + }; + + // for non-ghosts compute the derivative numerically + if (!this->elementIsGhost()) + { + static const NumericEpsilon eps_{this->problem().paramGroup()}; + static const int numDiffMethod = getParamFromGroup(this->problem().paramGroup(), "Assembly.NumericDifferenceMethod"); + NumericDifferentiation::partialDerivative(evalSource, elemSol[0][pvIdx], partialDeriv, sourceResidual, + eps_(elemSol[0][pvIdx], pvIdx), numDiffMethod); + } + + // for ghost elements we assemble a 1.0 where the primary variable and zero everywhere else + // as we always solve for a delta of the solution with repect to the initial solution this + // results in a delta of zero for ghosts + else partialDeriv[pvIdx] = 1.0; + + // add the current partial derivatives to the global jacobian matrix + // only diagonal entries for the source term + for (int eqIdx = 0; eqIdx < numEq; eqIdx++) + A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx]; + + // restore the original state of the scv's volume variables + curVolVars = origVolVars; + + // restore the current element solution + elemSol[0][pvIdx] = origPriVars[pvIdx]; + + // restore the undeflected state of the coupling context + this->couplingManager().updateCouplingContext(domainI, *this, domainI, globalI, elemSol[0], pvIdx); + } + + // restore original state of the flux vars cache in case of global caching. + // This has to be done in order to guarantee that everything is in an undeflected + // state before the assembly of another element is called. In the case of local caching + // this is obsolete because the elemFluxVarsCache used here goes out of scope after this. + // We only have to do this for the last primary variable, for all others the flux var cache + // is updated with the correct element volume variables before residual evaluations + if (enableGridFluxVarsCache) + gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars); + } + + // analytic diff for source + else + this->localResidual().addSourceDerivatives(A[globalI][globalI], problem, element, fvGeometry, curElemVolVars, scv); + // add flux derivatives for each scvf for (const auto& scvf : scvfs(fvGeometry)) @@ -859,8 +936,14 @@ public: return residual; } else + { + // compute potential additional derivatives caused by the coupling + const auto origResidual = this->evalLocalResidual(); + this->couplingManager().evalAdditionalDomainDerivatives(domainI, *this, origResidual, A, gridVariables); + // return element residual - return this->evalLocalResidual()[0]; + return origResidual[0]; + } } /*! @@ -877,9 +960,12 @@ public: // Calculate derivatives of all dofs in the element // // with respect to all dofs in the coupling stencil. // /////////////////////////////////////////////////////// - this->couplingManager().addCouplingDerivatives(A, domainI, *this, domainJ); + static const bool useNumDiff = getParam("Assembly.UseNumericDifferentiationForCouplingDerivatives", false); + if (useNumDiff) + assembleJacobianCouplingNumeric(domainJ, A, res, gridVariables); + else + this->couplingManager().addCouplingDerivatives(A, domainI, *this, domainJ); - // add the current partial derivatives to the global jacobian matrix if constexpr (Problem::enableInternalDirichletConstraints()) { @@ -896,6 +982,106 @@ public: } } + /*! + * \brief Computes the derivatives with respect to the given element and adds them + * to the global matrix. + * \todo This is just a copy of the assembleJacobianCoupling of the numeric-diff local assembler + */ + template + void assembleJacobianCouplingNumeric(Dune::index_constant domainJ, JacobianBlock& A, + const LocalResidualValues& res, GridVariables& gridVariables) + { + //////////////////////////////////////////////////////////////////////////////////////////////////////// + // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. // + //////////////////////////////////////////////////////////////////////////////////////////////////////// + + // get some aliases for convenience + const auto& element = this->element(); + const auto& fvGeometry = this->fvGeometry(); + const auto& fvGridGeometry = fvGeometry.fvGridGeometry(); + auto&& curElemVolVars = this->curElemVolVars(); + auto&& elemFluxVarsCache = this->elemFluxVarsCache(); + + // get stencil informations + const auto globalI = fvGridGeometry.elementMapper().index(element); + const auto& stencil = this->couplingManager().couplingStencil(domainI, element, domainJ); + const auto& curSolJ = this->curSol()[domainJ]; + + // make sure the flux variables cache does not contain any artifacts from prior differentiation + elemFluxVarsCache.update(element, fvGeometry, curElemVolVars); + + // convenience lambda for call to update self + auto updateCoupledVariables = [&] () + { + // Update ourself after the context has been modified. Depending on the + // type of caching, other objects might have to be updated. All ifs can be optimized away. + if (enableGridFluxVarsCache) + { + if (enableGridVolVarsCache) + { + this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), gridVariables.gridFluxVarsCache()); + this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), elemFluxVarsCache); + } + else + { + this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, gridVariables.gridFluxVarsCache()); + this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, elemFluxVarsCache); + } + } + else + { + if (enableGridVolVarsCache) + this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), elemFluxVarsCache); + else + this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, elemFluxVarsCache); + } + }; + + for (const auto& globalJ : stencil) + { + // undeflected privars and privars to be deflected + const auto origPriVarsJ = curSolJ[globalJ]; + auto priVarsJ = origPriVarsJ; + + const auto origResidual = this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ)[0]; + + for (int pvIdx = 0; pvIdx < JacobianBlock::block_type::cols; ++pvIdx) + { + auto evalCouplingResidual = [&](Scalar priVar) + { + // update the volume variables and the flux var cache + priVarsJ[pvIdx] = priVar; + this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx); + updateCoupledVariables(); + return this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ)[0]; + }; + + // derive the residuals numerically + LocalResidualValues partialDeriv(0.0); + const auto& paramGroup = this->assembler().problem(domainJ).paramGroup(); + static const int numDiffMethod = getParamFromGroup(paramGroup, "Assembly.NumericDifferenceMethod"); + static const auto epsCoupl = this->couplingManager().numericEpsilon(domainJ, paramGroup); + NumericDifferentiation::partialDerivative(evalCouplingResidual, priVarsJ[pvIdx], partialDeriv, origResidual, + epsCoupl(priVarsJ[pvIdx], pvIdx), numDiffMethod); + + // add the current partial derivatives to the global jacobian matrix + for (int eqIdx = 0; eqIdx < numEq; eqIdx++) + A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx]; + + // restore the current element solution + priVarsJ[pvIdx] = origPriVarsJ[pvIdx]; + + // restore the undeflected state of the coupling context + this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx); + } + + // restore original state of the flux vars cache and/or vol vars in case of global caching. + // This has to be done in order to guarantee that the undeflected residual computation done + // above is correct when jumping to the next coupled element. + updateCoupledVariables(); + } + } + }; // implicit assebler with analytic Jacobian } // end namespace Dumux -- GitLab From 839e2ad52bca9723e75ff8239badc4cc3d97ef4d Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Thu, 29 Aug 2019 21:17:35 +0200 Subject: [PATCH 07/10] [test][multidomain] Add richards_1p test with parital analytic diff --- .../embedded/1d3d/1p_richards/CMakeLists.txt | 72 ++++++++++++++++++- .../embedded/1d3d/1p_richards/main.cc | 19 +++-- .../embedded/1d3d/1p_richards/params.input | 4 ++ .../embedded/1d3d/1p_richards/problem_root.hh | 5 +- 4 files changed, 91 insertions(+), 9 deletions(-) diff --git a/test/multidomain/embedded/1d3d/1p_richards/CMakeLists.txt b/test/multidomain/embedded/1d3d/1p_richards/CMakeLists.txt index 690280bb66..fadd9441b3 100644 --- a/test/multidomain/embedded/1d3d/1p_richards/CMakeLists.txt +++ b/test/multidomain/embedded/1d3d/1p_richards/CMakeLists.txt @@ -2,6 +2,7 @@ dumux_add_test(NAME test_md_embedded_1d3d_1p_richards_tpfatpfa LABELS multidomain multidomain_embedded 1p richards SOURCES main.cc COMPILE_DEFINITIONS SOILTYPETAG=SoilCC + COMPILE_DEFINITIONS DIFFMETHOD=DiffMethod::numeric COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py CMAKE_GUARD dune-foamgrid_FOUND TIMEOUT 1500 @@ -13,10 +14,27 @@ dumux_add_test(NAME test_md_embedded_1d3d_1p_richards_tpfatpfa --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfatpfa params.input \ -Vtk.OutputName test_md_embedded_1d3d_1p_richards_tpfatpfa") -dumux_add_test(NAME test_md_embedded_1d3d_1p_richards_tpfabox +dune_add_test(NAME test_md_embedded_1d3d_1p_richards_tpfatpfa_anadiff + LABELS multidomain multidomain_embedded 1p richards + SOURCES main.cc + COMPILE_DEFINITIONS SOILTYPETAG=SoilCC + COMPILE_DEFINITIONS DIFFMETHOD=DiffMethod::analytic + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMAKE_GUARD dune-foamgrid_FOUND + TIMEOUT 1500 + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_tpfatpfa_1d-reference.vtp + ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfatpfa_anadiff_1d-00004.vtp + ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_tpfatpfa_3d-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfatpfa_anadiff_3d-00004.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfatpfa_anadiff params.input \ + -Vtk.OutputName test_md_embedded_1d3d_1p_richards_tpfatpfa_anadiff") + +dune_add_test(NAME test_md_embedded_1d3d_1p_richards_tpfabox LABELS multidomain multidomain_embedded 1p richards SOURCES main.cc COMPILE_DEFINITIONS SOILTYPETAG=SoilBox + COMPILE_DEFINITIONS DIFFMETHOD=DiffMethod::numeric COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py CMAKE_GUARD dune-foamgrid_FOUND TIMEOUT 1500 @@ -26,5 +44,55 @@ dumux_add_test(NAME test_md_embedded_1d3d_1p_richards_tpfabox ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_tpfabox_3d-reference.vtu ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfabox_3d-00004.vtu --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfabox params.input \ - -Vtk.OutputName test_md_embedded_1d3d_1p_richards_tpfabox") + -Vtk.OutputName test_md_embedded_1d3d_1p_richards_tpfabox") + +dune_add_test(NAME test_md_embedded_1d3d_1p_richards_tpfabox_anadiff + LABELS multidomain multidomain_embedded 1p richards + SOURCES main.cc + COMPILE_DEFINITIONS SOILTYPETAG=SoilBox + COMPILE_DEFINITIONS DIFFMETHOD=DiffMethod::analytic + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMAKE_GUARD dune-foamgrid_FOUND + TIMEOUT 1500 + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_tpfatpfa_1d-reference.vtp + ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfabox_anadiff_1d-00004.vtp + ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_tpfabox_3d-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfabox_anadiff_3d-00004.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfabox_anadiff params.input \ + -Vtk.OutputName test_md_embedded_1d3d_1p_richards_tpfabox_anadiff") + +dune_add_test(NAME test_md_embedded_1d3d_1p_richards_boxbox + LABELS multidomain multidomain_embedded 1p richards + SOURCES main.cc + COMPILE_DEFINITIONS SOILTYPETAG=SoilBox + COMPILE_DEFINITIONS ROOTTYPETAG=RootBox + COMPILE_DEFINITIONS DIFFMETHOD=DiffMethod::numeric + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMAKE_GUARD dune-foamgrid_FOUND + TIMEOUT 1500 + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_tpfatpfa_1d-reference.vtp + ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_boxbox_1d-00004.vtp + ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_boxbox_3d-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_boxbox_3d-00004.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_boxbox params.input \ + -Vtk.OutputName test_md_embedded_1d3d_1p_richards_boxbox") + +dune_add_test(NAME test_md_embedded_1d3d_1p_richards_boxbox_anadiff + LABELS multidomain multidomain_embedded 1p richards + SOURCES main.cc + COMPILE_DEFINITIONS SOILTYPETAG=SoilBox + COMPILE_DEFINITIONS ROOTTYPETAG=RootBox + COMPILE_DEFINITIONS DIFFMETHOD=DiffMethod::analytic + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMAKE_GUARD dune-foamgrid_FOUND + TIMEOUT 1500 + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_tpfatpfa_1d-reference.vtp + ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_boxbox_anadiff_1d-00004.vtp + ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_boxbox_3d-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_boxbox_anadiff_3d-00004.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_boxbox_anadiff params.input \ + -Vtk.OutputName test_md_embedded_1d3d_1p_richards_boxbox_anadiff") dune_symlink_to_source_files(FILES "params.input") diff --git a/test/multidomain/embedded/1d3d/1p_richards/main.cc b/test/multidomain/embedded/1d3d/1p_richards/main.cc index dce774e341..4a99049f63 100644 --- a/test/multidomain/embedded/1d3d/1p_richards/main.cc +++ b/test/multidomain/embedded/1d3d/1p_richards/main.cc @@ -25,6 +25,13 @@ #include +#ifndef SOILTYPETAG +#define SOILTYPETAG SoilCC +#endif +#ifndef ROOTTYPETAG +#define ROOTTYPETAG RootCC +#endif + #include #include #include @@ -56,12 +63,12 @@ namespace Properties { template struct CouplingManager { - using Traits = MultiDomainTraits; + using Traits = MultiDomainTraits; using type = EmbeddedCouplingManager1d3d; }; template -struct CouplingManager +struct CouplingManager { using Traits = MultiDomainTraits; using type = EmbeddedCouplingManager1d3d; @@ -70,11 +77,11 @@ struct CouplingManager template struct PointSource { using type = typename GetPropType::PointSourceTraits::template PointSource<0>; }; template -struct PointSource { using type = typename GetPropType::PointSourceTraits::template PointSource<1>; }; +struct PointSource { using type = typename GetPropType::PointSourceTraits::template PointSource<1>; }; template struct PointSourceHelper { using type = typename GetPropType::PointSourceTraits::template PointSourceHelper<0>; }; template -struct PointSourceHelper { using type = typename GetPropType::PointSourceTraits::template PointSourceHelper<1>; }; +struct PointSourceHelper { using type = typename GetPropType::PointSourceTraits::template PointSourceHelper<1>; }; } // end namespace Properties } // end namespace Dumux @@ -95,7 +102,7 @@ int main(int argc, char** argv) // Define the sub problem type tags using BulkTypeTag = Properties::TTag::SOILTYPETAG; - using LowDimTypeTag = Properties::TTag::Root; + using LowDimTypeTag = Properties::TTag::ROOTTYPETAG; // try to create a grid (from the given grid file or the input file) // for both sub-domains @@ -188,7 +195,7 @@ int main(int argc, char** argv) timeLoop->setMaxTimeStepSize(maxDt); // the assembler with time loop for instationary problem - using Assembler = MultiDomainFVAssembler; + using Assembler = MultiDomainFVAssembler; auto assembler = std::make_shared(std::make_tuple(bulkProblem, lowDimProblem), std::make_tuple(bulkFvGridGeometry, lowDimFvGridGeometry), std::make_tuple(bulkGridVariables, lowDimGridVariables), diff --git a/test/multidomain/embedded/1d3d/1p_richards/params.input b/test/multidomain/embedded/1d3d/1p_richards/params.input index fe09ff1a66..408a900075 100644 --- a/test/multidomain/embedded/1d3d/1p_richards/params.input +++ b/test/multidomain/embedded/1d3d/1p_richards/params.input @@ -35,3 +35,7 @@ Problem.Name = 1d InitialSoilPressure = -0.9429e4 # [Pa] InitialRootPressure = -1.2e6 # [Pa] TranspirationRate = 2.15e-8 # [kg / s] + +[Assembly] +UseNumericDifferentiationForSourceDerivatives = true +UseNumericDifferentiationForCouplingDerivatives = true diff --git a/test/multidomain/embedded/1d3d/1p_richards/problem_root.hh b/test/multidomain/embedded/1d3d/1p_richards/problem_root.hh index fd2cf05928..801137810c 100644 --- a/test/multidomain/embedded/1d3d/1p_richards/problem_root.hh +++ b/test/multidomain/embedded/1d3d/1p_richards/problem_root.hh @@ -32,6 +32,7 @@ #include #include #include +#include #include #include @@ -50,7 +51,9 @@ namespace Properties { // Create new type tags namespace TTag { -struct Root { using InheritsFrom = std::tuple; }; +struct Root { using InheritsFrom = std::tuple; }; +struct RootCC { using InheritsFrom = std::tuple; }; +struct RootBox { using InheritsFrom = std::tuple; }; } // end namespace TTag // Set the grid type -- GitLab From 247f57b09e6f026991e955cdf1a87772c143645c Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 2 Oct 2020 10:41:55 +0200 Subject: [PATCH 08/10] Fixup problem interface --- dumux/common/fvproblem.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh index 6b27cc5555..49118f70bf 100644 --- a/dumux/common/fvproblem.hh +++ b/dumux/common/fvproblem.hh @@ -415,7 +415,7 @@ public: * \brief Add source term derivative to the Jacobian * \note Only needed in case of analytic differentiation and solution dependent sources */ - template + template void addSourceDerivatives(MatrixBlock& block, const Element& element, const FVElementGeometry& fvGeometry, -- GitLab From fd0895fb22bbc747d817a23919b83547e954a733 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 2 Oct 2020 10:42:24 +0200 Subject: [PATCH 09/10] ficup grid geoetmry --- dumux/multidomain/subdomainboxlocalassembler.hh | 2 +- dumux/multidomain/subdomaincclocalassembler.hh | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/dumux/multidomain/subdomainboxlocalassembler.hh b/dumux/multidomain/subdomainboxlocalassembler.hh index df9e7bbbac..3b599e8e03 100644 --- a/dumux/multidomain/subdomainboxlocalassembler.hh +++ b/dumux/multidomain/subdomainboxlocalassembler.hh @@ -713,7 +713,7 @@ class SubDomainBoxLocalAssembler::numEq(); static constexpr bool enableGridFluxVarsCache = getPropValue(); static constexpr bool enableGridVolVarsCache = getPropValue(); - enum { dim = GetPropType::dimension }; + enum { dim = GetPropType::GridView::dimension }; public: using ParentType::ParentType; diff --git a/dumux/multidomain/subdomaincclocalassembler.hh b/dumux/multidomain/subdomaincclocalassembler.hh index 501dc85639..3c455b868a 100644 --- a/dumux/multidomain/subdomaincclocalassembler.hh +++ b/dumux/multidomain/subdomaincclocalassembler.hh @@ -795,12 +795,12 @@ public: const auto& problem = this->problem(); const auto& element = this->element(); const auto& fvGeometry = this->fvGeometry(); - const auto& fvGridGeometry = fvGeometry.fvGridGeometry(); + const auto& gridGeometry = fvGeometry.gridGeometry(); auto&& curElemVolVars = this->curElemVolVars(); const auto& elemFluxVarsCache = this->elemFluxVarsCache(); // get reference to the element's current vol vars - const auto globalI = fvGridGeometry.elementMapper().index(element); + const auto globalI = gridGeometry.elementMapper().index(element); const auto& scv = fvGeometry.scv(globalI); const auto& volVars = curElemVolVars[scv]; @@ -822,7 +822,7 @@ public: const auto origVolVars = curVolVars; // element solution container to be deflected - auto elemSol = elementSolution(element, curSol, fvGridGeometry); + auto elemSol = elementSolution(element, curSol, gridGeometry); // derivatives in the neighbors with repect to the current elements LocalResidualValues partialDeriv; @@ -998,12 +998,12 @@ public: // get some aliases for convenience const auto& element = this->element(); const auto& fvGeometry = this->fvGeometry(); - const auto& fvGridGeometry = fvGeometry.fvGridGeometry(); + const auto& gridGeometry = fvGeometry.gridGeometry(); auto&& curElemVolVars = this->curElemVolVars(); auto&& elemFluxVarsCache = this->elemFluxVarsCache(); // get stencil informations - const auto globalI = fvGridGeometry.elementMapper().index(element); + const auto globalI = gridGeometry.elementMapper().index(element); const auto& stencil = this->couplingManager().couplingStencil(domainI, element, domainJ); const auto& curSolJ = this->curSol()[domainJ]; -- GitLab From 0310c901efd91fbd0ccb43bf94ab8b5fe3f7f980 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 2 Oct 2020 10:42:50 +0200 Subject: [PATCH 10/10] ficup dirichlet contraints --- dumux/multidomain/subdomaincclocalassembler.hh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/dumux/multidomain/subdomaincclocalassembler.hh b/dumux/multidomain/subdomaincclocalassembler.hh index 3c455b868a..3882f22361 100644 --- a/dumux/multidomain/subdomaincclocalassembler.hh +++ b/dumux/multidomain/subdomaincclocalassembler.hh @@ -775,7 +775,6 @@ class SubDomainCCLocalAssembler(); - static constexpr int numEq = GetPropType::numEq(); static constexpr bool enableGridFluxVarsCache = getPropValue(); static constexpr bool enableGridVolVarsCache = getPropValue(); @@ -968,10 +967,11 @@ public: if constexpr (Problem::enableInternalDirichletConstraints()) { - + const auto globalI = this->fvGeometry().gridGeometry().elementMapper().index(this->element()); + const auto& stencil = this->couplingManager().couplingStencil(domainI, this->element(), domainJ); for (const auto globalJ : stencil) { - const auto& scv = fvGeometry.scv(globalI); + const auto& scv = this->fvGeometry().scv(globalI); const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv); if (internalDirichletConstraints.any()) for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) -- GitLab