From fe7a21a4213bd1fec334c5e2cccd733240384f69 Mon Sep 17 00:00:00 2001 From: Lukas Riedel Date: Fri, 16 Feb 2018 19:14:58 +0100 Subject: [PATCH] Finalize DG scheme options in local operator - fixed full upwdining: only use satCond for weighting - re-implemented option to disable weights in skeleton term - moved experimental settings function to Operator namespace --- dune/dorie/interface/simulation.cc | 13 +++- dune/dorie/interface/simulation.hh | 40 ---------- dune/dorie/solver/operator_DG.hh | 116 +++++++++++++++++++++-------- 3 files changed, 93 insertions(+), 76 deletions(-) diff --git a/dune/dorie/interface/simulation.cc b/dune/dorie/interface/simulation.cc index f98dad44..b764481d 100644 --- a/dune/dorie/interface/simulation.cc +++ b/dune/dorie/interface/simulation.cc @@ -32,9 +32,15 @@ Simulation::Simulation (Dune::MPIHelper& _helper, std::shared_ptr // --- Local Operators --- #ifdef EXPERIMENTAL_DG_FEATURES - const auto settings = read_experimental_operator_settings(inifile); + // read experimental settings from inifile + namespace OP = Dune::Dorie::Operator; + const auto settings = OP::read_experimental_operator_settings(inifile); + const auto method = std::get(settings); + const auto upwinding = std::get(settings); + const auto weights = std::get(settings); + slop = std::make_unique(inifile, *param, gv, *fboundary, *fsource, - settings.first, settings.second); + method, upwinding, weights); #else slop = std::make_unique(inifile,*param,gv,*fboundary,*fsource); #endif // EXPERIMENTAL_DG_FEATURES @@ -135,6 +141,7 @@ bool Simulation::compute_time_step () const RF t = controller->getTime(); const RF dt = controller->getDT(); bool exception = false; + try { // solve in parallel @@ -147,7 +154,7 @@ bool Simulation::compute_time_step () // solve sequentially else { pdesolver_seq->setMaxIterations(controller->getIterations()); - osm_seq->apply(t,dt,*uold,*unew); + osm_seq->apply(t, dt, *uold, *unew); } *uold = *unew; diff --git a/dune/dorie/interface/simulation.hh b/dune/dorie/interface/simulation.hh index 11b6f313..6683164d 100644 --- a/dune/dorie/interface/simulation.hh +++ b/dune/dorie/interface/simulation.hh @@ -127,46 +127,6 @@ protected: const int verbose; -private: - - /// Read experimental operator settings from a parameter file - /* \return Pair of RichardsDGMethod and RichardsDGUpwinding - * to be inserted into the local operator constructor. - */ - auto read_experimental_operator_settings (const Dune::ParameterTree& inifile) - -> std::pair - { - using namespace Operator; - RichardsDGMethod::Type method; - RichardsDGUpwinding::Type upwinding; - - const auto method_str = inifile.get("dg.method"); - if (method_str == "SIPG") - method = RichardsDGMethod::SIPG; - else if (method_str == "NIPG") - method = RichardsDGMethod::NIPG; - else if (method_str == "OOB") - method = RichardsDGMethod::OOB; - else if (method_str == "IIP") - method = RichardsDGMethod::IIP; - else - DUNE_THROW(Dune::IOError, "Unknown DG method '" + method_str + "'!"); - - const auto upwinding_str = inifile.get("dg.upwinding"); - if (upwinding_str == "none") - upwinding = RichardsDGUpwinding::none; - else if (upwinding_str == "semiUpwind") - upwinding = RichardsDGUpwinding::semiUpwind; - else if (upwinding_str == "fullUpwind") - upwinding = RichardsDGUpwinding::fullUpwind; - else - DUNE_THROW(Dune::IOError, "Unknown upwdinding '" + upwinding_str + "'!"); - - return std::make_pair(method, upwinding); - } - - public: /// Execute the simulation until tEnd is reached. diff --git a/dune/dorie/solver/operator_DG.hh b/dune/dorie/solver/operator_DG.hh index 3b5dce68..a857e13a 100755 --- a/dune/dorie/solver/operator_DG.hh +++ b/dune/dorie/solver/operator_DG.hh @@ -35,6 +35,15 @@ struct RichardsDGUpwinding }; }; +//! Gradient weighting +struct RichardsDGWeights +{ + enum Type { + weightsOn, //!< harmonic weighting of flux contributions + weightsOff //!< no weighting of flux contributions + }; +}; + //! Operator internal BC handling struct BCType { @@ -44,6 +53,52 @@ struct BCType }; }; +/// Read experimental operator settings from a parameter file +/* \return Tuple of RichardsDGMethod, RichardsDGUpwinding, and + * RichardsDGWeights to be inserted into the local operator constructor. + */ +inline auto read_experimental_operator_settings( + const Dune::ParameterTree &inifile) + -> std::tuple +{ + RichardsDGMethod::Type method; + RichardsDGUpwinding::Type upwinding; + RichardsDGWeights::Type weights; + + const auto method_str = inifile.get("dg.experimental.method"); + if (method_str == "SIPG") + method = RichardsDGMethod::SIPG; + else if (method_str == "NIPG") + method = RichardsDGMethod::NIPG; + else if (method_str == "OOB") + method = RichardsDGMethod::OOB; + else if (method_str == "IIP") + method = RichardsDGMethod::IIP; + else + DUNE_THROW(Dune::IOError, "Unknown DG method '" + method_str + "'!"); + + const auto upwinding_str + = inifile.get("dg.experimental.upwinding"); + if (upwinding_str == "none") + upwinding = RichardsDGUpwinding::none; + else if (upwinding_str == "semiUpwind") + upwinding = RichardsDGUpwinding::semiUpwind; + else if (upwinding_str == "fullUpwind") + upwinding = RichardsDGUpwinding::fullUpwind; + else + DUNE_THROW(Dune::IOError, "Unknown upwinding '" + upwinding_str + "'!"); + + const auto weights_str = inifile.get("dg.experimental.weights"); + if (weights_str) + weights = RichardsDGWeights::weightsOn; + else + weights = RichardsDGWeights::weightsOff; + + return std::make_tuple(method, upwinding, weights); +} + /// A local operator for solving the Richards equation using upwinding. /** Written by Ole Klein, IWR, Heidelberg and modified by the * DORiE developers, IUP, Heidelberg \n @@ -87,8 +142,9 @@ protected: const Parameter& param; //!< class providing equation parameters const Boundary& boundary; //!< class providing boundary conditions const SourceTerm& sourceTerm; //!< class providing source term information - RichardsDGMethod::Type method; //!< Interior penalty type - RichardsDGUpwinding::Type upwinding; //!< DG method weights switch + const RichardsDGMethod::Type method; //!< Interior penalty type + const RichardsDGUpwinding::Type upwinding; //!< DG method weights switch + const RichardsDGWeights::Type weights; //!< Gradient weighting const Mapper mapper; //!< maps grid elements (cells) to natural numbers const int intorderadd; //!< integration order addend @@ -120,6 +176,7 @@ public: * \param boundary_ Object prividing boundary conditions * \param sourceTerm_ Object providing source term information * \param method_ DG method type \see RichardsDGMethod::Type + * \param upwinding_ Upwinding type \see RichardsDGUpwinding::Type * \param weights_ DG weigths switch \see RichardsDGWeights::Type * \param intorderadd_ Addend to integration order * \param quadrature_factor_ Factor to FEM integration order @@ -128,11 +185,12 @@ public: Boundary& boundary_, const SourceTerm& sourceTerm_, RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG, RichardsDGUpwinding::Type upwinding_ = RichardsDGUpwinding::semiUpwind, + RichardsDGWeights::Type weights_ = RichardsDGWeights::weightsOn, int intorderadd_ = 2, int quadrature_factor_ = 2) : Dune::PDELab::NumericalJacobianVolume >(1.e-7), Dune::PDELab::NumericalJacobianSkeleton >(1.e-7), Dune::PDELab::NumericalJacobianBoundary >(1.e-7), - param(param_), boundary(boundary_), sourceTerm(sourceTerm_), method(method_), upwinding(upwinding_), mapper(view_), + param(param_), boundary(boundary_), sourceTerm(sourceTerm_), method(method_), upwinding(upwinding_), weights(weights_), mapper(view_), intorderadd(intorderadd_), quadrature_factor(quadrature_factor_), penalty_factor(config.get("dg.penaltyFactor")), time(0.0), @@ -330,18 +388,21 @@ public: // compute weights RF omega_s = 0.5, omega_n = 0.5; - if (upwinding == RichardsDGUpwinding::none - || upwinding == RichardsDGUpwinding::fullUpwind) - { - omega_s = relCond_n * satCond_n - / ((relCond_s * satCond_s) + (relCond_n * satCond_n)); - omega_n = relCond_s * satCond_s - / ((relCond_s * satCond_s) + (relCond_n * satCond_n)); - } - else if (upwinding == RichardsDGUpwinding::semiUpwind) + if (weights == RichardsDGWeights::weightsOn) { - omega_s = satCond_n / (satCond_s + satCond_n); - omega_n = satCond_s / (satCond_s + satCond_n); + if (upwinding == RichardsDGUpwinding::none) + { + omega_s = relCond_n * satCond_n + / ((relCond_s * satCond_s) + (relCond_n * satCond_n)); + omega_n = relCond_s * satCond_s + / ((relCond_s * satCond_s) + (relCond_n * satCond_n)); + } + else if (upwinding == RichardsDGUpwinding::semiUpwind + || upwinding == RichardsDGUpwinding::fullUpwind) + { + omega_s = satCond_n / (satCond_s + satCond_n); + omega_n = satCond_s / (satCond_s + satCond_n); + } } // integration factor @@ -352,16 +413,15 @@ public: // compute numerical flux RF numFlux_s = satCond_s * ( - (gradu_s * normal) + penalty * jump); RF numFlux_n = satCond_n * ( - (gradu_n * normal) + penalty * jump); - if (upwinding == RichardsDGUpwinding::none - || upwinding == RichardsDGUpwinding::fullUpwind) + if (upwinding == RichardsDGUpwinding::none) { numFlux_s *= relCond_s; numFlux_n *= relCond_n; } const RF numFlux = omega_s * numFlux_s + omega_n * numFlux_n; - // compute upwind conductivity factor - RF relCond_upwind = 1.0; + // compute upwind quantities + RF relCond_upwind, satCond_upwind, numFlux_upwind; if (upwinding == RichardsDGUpwinding::semiUpwind || upwinding == RichardsDGUpwinding::fullUpwind) { @@ -370,26 +430,16 @@ public: const RF head_upwind = param.headMillerToRef(u_s, p_global_s); const RF saturation_upwind = param.saturation(head_upwind, p_global_s); relCond_upwind = param.condFactor(saturation_upwind, p_global_s); + satCond_upwind = satCond_s; + numFlux_upwind = numFlux_s; } else { const RF head_upwind = param.headMillerToRef(u_n, p_global_n); const RF saturation_upwind = param.saturation(head_upwind, p_global_n); relCond_upwind = param.condFactor(saturation_upwind, p_global_n); - } - } - - // compute upwind flux - RF numFlux_upwind; - RF satCond_upwind; - if (upwinding == RichardsDGUpwinding::fullUpwind) { - if (numFlux > 0) { - numFlux_upwind = numFlux_s; - satCond_upwind = satCond_s; - } - else { - numFlux_upwind = numFlux_n; satCond_upwind = satCond_n; + numFlux_upwind = numFlux_n; } } @@ -414,9 +464,9 @@ public: else if (upwinding == RichardsDGUpwinding::fullUpwind) { for (unsigned int i = 0; i < lfsv_s.size(); i++) - r_s.accumulate(lfsv_s, i, numFlux_upwind * phi_s[i] * factor); + r_s.accumulate(lfsv_s, i, relCond_upwind * numFlux_upwind * phi_s[i] * factor); for (unsigned int i = 0; i < lfsv_n.size(); i++) - r_n.accumulate(lfsv_n, i, - numFlux_upwind * phi_n[i] * factor); + r_n.accumulate(lfsv_n, i, - relCond_upwind * numFlux_upwind * phi_n[i] * factor); } // update residual (symmetry term) -- GitLab