From 8ce31d123192c95f5225626ffc5e2e0f98593589 Mon Sep 17 00:00:00 2001 From: Santiago Ospina Date: Tue, 11 Sep 2018 18:18:44 +0200 Subject: [PATCH] Added a general flux reconstruction method, unfortunately Raviart Thomas elements are not implemented in all of the dimensions and degrees we need Signed-off-by: Santiago Ospina --- dune/dorie/interface/richards_simulation.cc | 15 ++ dune/dorie/solver/util_flux_reconstruction.hh | 158 ++++++++++++++++++ dune/dorie/solver/util_writer.hh | 5 +- python/dorie/dorie/cli/cmds.py | 4 +- 4 files changed, 178 insertions(+), 4 deletions(-) create mode 100644 dune/dorie/solver/util_flux_reconstruction.hh diff --git a/dune/dorie/interface/richards_simulation.cc b/dune/dorie/interface/richards_simulation.cc index f8a503f5..e2518288 100644 --- a/dune/dorie/interface/richards_simulation.cc +++ b/dune/dorie/interface/richards_simulation.cc @@ -1,4 +1,5 @@ #include "richards_simulation.hh" +#include "../solver/util_flux_reconstruction.hh" namespace Dune{ namespace Dorie{ @@ -213,18 +214,32 @@ void RichardsSimulation::write_data (RF time) if constexpr (Traits::write_output) { update_adapters(); + + std::shared_ptr> lift = std::make_shared< + Dune::Dorie::FluxReconstruction>(udgf); + if (inifile.get("output.vertexData")) { vtkwriter->addVertexData(udgf,"head"); vtkwriter->addVertexData(fluxdgf,"flux"); vtkwriter->addVertexData(condgf,"K_0"); vtkwriter->addVertexData(waterdgf,"theta_w"); vtkwriter->addVertexData(satdgf,"Theta"); + // vtkwriter->addVertexData(lift,"lift"); } else { vtkwriter->addCellData(udgf,"head"); vtkwriter->addCellData(fluxdgf,"flux"); vtkwriter->addCellData(condgf,"K_0"); vtkwriter->addCellData(waterdgf,"theta_w"); vtkwriter->addCellData(satdgf,"Theta"); + // vtkwriter->addCellData(lift,"lift"); } try{ vtkwriter->write(time,output_type); diff --git a/dune/dorie/solver/util_flux_reconstruction.hh b/dune/dorie/solver/util_flux_reconstruction.hh new file mode 100644 index 00000000..037248c7 --- /dev/null +++ b/dune/dorie/solver/util_flux_reconstruction.hh @@ -0,0 +1,158 @@ +// -*- tab-width: 4; indent-tabs-mode: nil -*- +#ifndef DUNE_DORIE_FLUX_RECONSTRUCTION_HH +#define DUNE_DORIE_FLUX_RECONSTRUCTION_HH + +#include + +namespace Dune { +namespace Dorie { + +/** + * @brief Class for global lifting. + * + * @tparam { description } + */ +template +class GlobalLifting + : public Dune::PDELab::LocalOperatorDefaultFlags +{ +public: + GlobalLifting(std::shared_ptr grid_function) + : _grid_function(grid_function) + {} + + enum {doAlphaSkeleton = true}; + enum {doSkeletonTwoSided = true}; // FIXME, my programer was too lazy to do the two sides at once + + /*-----------------------------------------------------------------------*//** + * @brief Skeleton integral depending on test and ansatz functions. Each + * face is only visited once since this method is symmetric. + * + * @param[in] ig The intersection entity of the grid (inside + outside + * entities) + * @param[in] lfsu_i The inside ansatz local function space + * @param[in] x_i The coefficients of the lfsu_i + * @param[in] lfsv_i The inside test local function space + * @param[in] lfsu_o The outside ansatz local function space + * @param[in] x_o The coefficients of the lfsu_o + * @param[in] lfsv_o The outside test local function space + * @param r_i The view of the residual vector w.r.t lfsu_i + * @param r_o The view of the residual vector w.r.t lfsu_o + * + * @tparam IG The type for ig + * @tparam LFSU The type for lfsu_i and lfsu_o + * @tparam X The type for x_i and x_o + * @tparam LFSV The type for lfsv_i and lfsv_o + * @tparam R The type for r_i and r_o + */ + template + void alpha_skeleton(const IG& ig, + const LFSU& lfsu_i, const X& x_i, const LFSV& lfsv_i, + const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o, + R& r_i, R& r_o) const + { + + using RF = typename LFSU::Traits::FiniteElementType:: + Traits::LocalBasisType::Traits::RangeFieldType; + + const int intorder = 1; //FIXME + auto face_geo = ig.geometry(); + + assert(lfsu_i.size()==lfsv_i.size()); + assert(lfsu_o.size()==lfsv_o.size()); + + int n = lfsu_i.size(); + Dune::DynamicMatrix A(n,n); + Dune::DynamicVector b(n); + + // loop over quadrature points and integrate normal flux + for (const auto& it : quadratureRule(face_geo,intorder)) { + for (std::size_t i = 0; i < lfsv_i.size(); ++i) { + for (std::size_t j = 0; j < lfsv_o.size(); ++j) { + if (lfsv_i.localIndex(i) == lfsv_o.localIndex(j)) { + const auto average = 0.5*(x_i(lfsv_i,i) + x_o(lfsv_o,j)); // FIXME: average of a vector field! + const auto normal = ig.unitOuterNormal(it.position()); + typename GF::Traits::RangeType flux; + typename GF::Traits::DomainType x_local = ig.geometryInInside().global(it.position()); + _grid_function->evaluate(ig.inside(),x_local,flux); + b[i] += average * flux * normal ; // FIXME eval theta = jump in flux + } + } + for (std::size_t j = 0; j < lfsu_i.size(); ++j) + A[j][i] = x_i(lfsu_i,j)*x_i(lfsv_i,i); // Check order of i and j access + } + } + Dune::DynamicVector x(n); + A.solve(b,x); + for (std::size_t i = 0; i < lfsu_i.size(); ++i) + r_i.accumulate(lfsu_i,i,x[i]); + } +private: + std::shared_ptr _grid_function; +}; + +template +class FluxReconstruction : public Dune::PDELab::GridFunctionBase>,FluxReconstruction> +{ +public: + using Traits = Dune::PDELab::GridFunctionTraits>; +private: + static constexpr int dim = d; + using DF = typename Traits::DomainFieldType; + + using FEM = Dune::PDELab::RaviartThomasLocalFiniteElementMap; + using CON = Dune::PDELab::NoConstraints; + using VBE = Dune::PDELab::ISTL::VectorBackend<>; + using GFS = Dune::PDELab::GridFunctionSpace; + using X = Dune::PDELab::Backend::Vector; + using DGF = Dune::PDELab::DiscreteGridFunction; + +public: + template + FluxReconstruction(std::shared_ptr grid_function) + : _gv(grid_function->getGridView()) + , _fem(_gv) + , _gfs(_gv,_fem) + , _r(_gfs,0.0) + { + using MBE = Dune::PDELab::ISTL::BCRSMatrixBackend<>; + MBE mbe(5); + + using LOP = GlobalLifting; + LOP lop(grid_function); + using GO = Dune::PDELab::GridOperator; + GO go(_gfs,_gfs,lop,mbe); + X x(_gfs,1.0); + using LS = Dune::PDELab::ISTLBackend_SEQ_BCGS_AMG_SSOR; + LS ls; + go.residual(x,_r); + _dfg = std::make_unique(_gfs,_r); + } + + void evaluate ( const typename Traits::ElementType& e, + const typename Traits::DomainType& x, + typename Traits::RangeType& y) const + { + // _dfg->evaluate(e,x,y); + } + + const typename Traits::GridViewType& getGridView() const + { + return _gv; + } + +private: + GV _gv; + FEM _fem; + GFS _gfs; + X _r; + std::unique_ptr _dfg; +}; + +} // namespace Dorie +} // namespace Dune + + +#endif // DUNE_DORIE_FLUX_RECONSTRUCTION_HH diff --git a/dune/dorie/solver/util_writer.hh b/dune/dorie/solver/util_writer.hh index 81ff610c..55213e9f 100644 --- a/dune/dorie/solver/util_writer.hh +++ b/dune/dorie/solver/util_writer.hh @@ -44,7 +44,8 @@ public: static_assert(std::is_same::value, "GridFunctionVTKSequenceWriter only works with only one GridView type."); - auto vtk_gf_ptr = Dune::PDELab::makeVTKGridFunctionAdapter(gf_ptr,name); + + auto vtk_gf_ptr = std::make_shared >(gf_ptr, name); Dune::VTKSequenceWriter::addVertexData(vtk_gf_ptr); } @@ -62,7 +63,7 @@ public: static_assert(std::is_same::value, "GridFunctionVTKSequenceWriter only works with only one GridView type."); - auto vtk_gf_ptr = Dune::PDELab::makeVTKGridFunctionAdapter(gf_ptr,name); + auto vtk_gf_ptr = std::make_shared >(gf_ptr, name); Dune::VTKSequenceWriter::addCellData(vtk_gf_ptr); } diff --git a/python/dorie/dorie/cli/cmds.py b/python/dorie/dorie/cli/cmds.py index bea6ebbf..4d645e1f 100644 --- a/python/dorie/dorie/cli/cmds.py +++ b/python/dorie/dorie/cli/cmds.py @@ -9,8 +9,8 @@ except NameError: pass # paths set by cmake -DORIEDIR = "/Users/lriedel/dune/dorie/build-cmake" -MPIEXEC = "/usr/local/bin/mpiexec" +DORIEDIR = "/home/saospina/Codes/DUNE_INSTALL/Debug/dorie" +MPIEXEC = "/usr/bin/mpiexec" MPIEXEC_NUMPROC_FLAG = "-n" MPIEXEC_PREFLAG = "" MPIEXEC_POSTFLAGS = "" -- GitLab