[dune-pdelab] Finite Volume Method with slope reconstruction / alpha_volume_post_skeleton
Christian Engwer
christian.engwer at uni-muenster.de
Tue May 9 14:35:36 CEST 2017
Hi Gregor,
let me guess... you use the one-sided-assembly?! In this case the face
is only evaluated once and this is from the element with the lower
index (or was it the higher index?). In any case, this means that the
face might be skipped for your current cell and thus it is missing in
your limiter... If htis is the case you should switch to
two-sided-assembly.
Christian
On Tue, May 09, 2017 at 01:54:53PM +0200, Gregor Corbin wrote:
> Hi Christian,
>
> in my finite volume method with slope reconstruction, I implemented a local
> operator that computes the slopes.
>
> As we discussed before, it uses the alpha_skeleton method to compute slopes
> at each face, and stores this information. Then in the
> alpha_volume_post_skeleton method, it uses the previously computed slopes.
>
> But unfortunately, the alpha_volume_post_skeleton will not be executed after
> all the skeleton terms are done, but rather in between. To test this, I
> created a minimal 3x3 YASP-Grid example and let the methods of the local
> operator print their name and some info about the current cell.
>
> Is this behavior intended? If yes, it seems like I have to define two
> operators and apply them successively.
>
> I have attached the test program and its output. I use version 2.5 of all
> dune modules.
>
> Cheers,
>
> Gregor
>
>
> On 28.04.2017 20:29, Christian Engwer wrote:
> > Hi Gregor,
> >
> > > If I use a grid operator for this, how do I get the solution on neighboring
> > > cells? I suppose I could use the alpha_skeleton method to communicate this
> > > information and then doing the reconstruction in the
> > > alpha_volume_post_skeleton method. But this would mean a lot of storage and
> > > copying overhead.
> > I don't understand why you think you have to communicate information?!
> > You use alpha_skeleton to compute the slope w.r.t. each neighbor. You
> > can store this information in a stateful local operator and reuse the
> > information on the same cell in the post-skeleton.
> >
> > > Also, the reconstruction has to be done before each stage of the time
> > > stepping. So would I have to call the reconstruction operator from the
> > > preStage() method of my original local operator?
> > Yes.
> >
> > > Another question related to this: Are the intersections in a YASPGrid
> > > ordered in a consistent way? I.e. has intersections(grid_view,cell)[0]
> > > always the same normal?
> > Yes.
> >
> > Ciao
> > Christian
>
> Hello World! This is LocalOperatorTest.
> This is a sequential program.
> [INFO] Domain size = 3 3
> [INFO] #cells = 9
> alpha_boundary: cell inside 0 face 12
> alpha_boundary: cell inside 0 face 0
> alpha_volume_post_skeleton : cell = 0 center = -1 -1
> alpha_skeleton: cell inside 1 cell outside 0 face 13
> alpha_boundary: cell inside 1 face 1
> alpha_volume_post_skeleton : cell = 1 center = 0 -1
> alpha_skeleton: cell inside 2 cell outside 1 face 14
> alpha_boundary: cell inside 2 face 15
> alpha_boundary: cell inside 2 face 2
> alpha_volume_post_skeleton : cell = 2 center = 1 -1
> alpha_boundary: cell inside 3 face 16
> alpha_skeleton: cell inside 3 cell outside 0 face 3
> alpha_volume_post_skeleton : cell = 3 center = -1 0
> alpha_skeleton: cell inside 4 cell outside 3 face 17
> alpha_skeleton: cell inside 4 cell outside 1 face 4
> alpha_volume_post_skeleton : cell = 4 center = 0 0
> alpha_skeleton: cell inside 5 cell outside 4 face 18
> alpha_boundary: cell inside 5 face 19
> alpha_skeleton: cell inside 5 cell outside 2 face 5
> alpha_volume_post_skeleton : cell = 5 center = 1 0
> alpha_boundary: cell inside 6 face 20
> alpha_skeleton: cell inside 6 cell outside 3 face 6
> alpha_boundary: cell inside 6 face 9
> alpha_volume_post_skeleton : cell = 6 center = -1 1
> alpha_skeleton: cell inside 7 cell outside 6 face 21
> alpha_skeleton: cell inside 7 cell outside 4 face 7
> alpha_boundary: cell inside 7 face 10
> alpha_volume_post_skeleton : cell = 7 center = 0 1
> alpha_skeleton: cell inside 8 cell outside 7 face 22
> alpha_boundary: cell inside 8 face 23
> alpha_skeleton: cell inside 8 cell outside 5 face 8
> alpha_boundary: cell inside 8 face 11
> alpha_volume_post_skeleton : cell = 8 center = 1 1
> // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
> // vi: set et ts=4 sw=2 sts=2:
>
> #ifdef HAVE_CONFIG_H
> # include "config.h"
> #endif
>
> #include <iostream>
> #include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
> #include <dune/common/exceptions.hh> // We use exceptions
>
> #include <dune/grid/yaspgrid.hh>
>
> #include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
> #include <dune/grid/common/scsgmapper.hh>
>
> #include <dune/pdelab/finiteelementmap/opbfem.hh>
> #include <dune/pdelab/backend/istl.hh>
> #include <dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
> #include <dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
> #include <dune/pdelab/constraints/common/constraints.hh>
> /*
> * local operator for slope reconstruction
> * test version, that will only print out some info
> */
> template <typename GridView>
> class Slopes:
> public Dune::PDELab::NumericalJacobianApplySkeleton<Slopes<GridView> >,
> public Dune::PDELab::NumericalJacobianApplyBoundary<Slopes<GridView> >,
> public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<>,
> public Dune::PDELab::FullSkeletonPattern,
> public Dune::PDELab::LocalOperatorDefaultFlags
> {
>
> public:
> enum{doPatternVolume = false}; // No volume terms
> enum{doPatternSkeleton = false}; // Use Skeleton terms
>
> enum{doAlphaVolume = false};
> enum{doAlphaSkeleton = true};
> enum{doAlphaVolumePostSkeleton = true};
> enum{doAlphaBoundary = true};
>
> using FMapper =
> Dune::SingleCodimSingleGeomTypeMapper<GridView,1>;
> using CMapper =
> Dune::SingleCodimSingleGeomTypeMapper<GridView,0>;
>
> Slopes(const GridView &gridView):
> grid_view_(gridView),
> fmapper_(grid_view_),
> cmapper_(grid_view_)
> {}
>
> // reconstruct slopes here
> template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
> void alpha_volume_post_skeleton (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
> {
> // get cell
> const auto& cell = eg.entity();
>
> // cell center
> auto geo = eg.geometry();
>
>
> std::cout << "alpha_volume_post_skeleton : cell = " << cmapper_.index(cell) << " center = " << geo.center() << std::endl;
> }
>
> template <typename IG,
> typename LFSU,
> typename X,
> typename LFSV,
> typename R>
> void alpha_skeleton(const IG &ig,
> const LFSU &lfsu_s,const X &x_s,const LFSV &lfsv_s,
> const LFSU &lfsu_n,const X &x_n,const LFSV &lfsv_n,
> R &r_s,R &r_n)const
> {
> const auto& cell_inside = ig.inside();
> const auto& cell_outside = ig.outside();
>
> int idx = ig.indexInInside();
> std::cout << "alpha_skeleton:"
> <<" cell inside " << cmapper_.index(cell_inside)
> <<" cell outside " << cmapper_.index(cell_outside)
> <<" face " << fmapper_.subIndex(cell_inside,idx,1) << std::endl;
> }
>
> template <typename IG,
> typename LFSU,
> typename X,
> typename LFSV,
> typename R>
> void alpha_boundary(const IG &ig,
> const LFSU &lfsu_s,const X &x_s,const LFSV &lfsv_s,R &r_s)const
> {
> const auto& cell_inside = ig.inside();
>
> int idx = ig.indexInInside();
> std::cout << "alpha_boundary:"
> <<" cell inside " << cmapper_.index(cell_inside)
> <<" face " << fmapper_.subIndex(cell_inside,idx,1) << std::endl;
> }
>
> private:
>
> const GridView &grid_view_;
>
> FMapper fmapper_;
> CMapper cmapper_;
> };
>
>
>
>
> int main(int argc, char** argv)
> {
> try{
> // Maybe initialize MPI
> Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
> std::cout << "Hello World! This is LocalOperatorTest." << std::endl;
> if(Dune::MPIHelper::isFake)
> std::cout<< "This is a sequential program." << std::endl;
> else
> std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
> <<" processes!"<<std::endl;
>
> // build 3x3 yasp grid
> constexpr int dim = 2;
> using Grid = Dune::YaspGrid<dim,Dune::EquidistantOffsetCoordinates<double,dim>>;
> using GridView = typename Grid::LeafGridView;
>
> Dune::FieldVector<double,dim> lower_left_corner{{-1.5,-1.5}};
> Dune::FieldVector<double,dim> upper_right_corner{{1.5,1.5}};
> const std::array<int,dim> sizes{{3,3}};
> std::bitset<dim> periodic; // everything false
> unsigned int overlap = 1;
>
> Grid grid(lower_left_corner,upper_right_corner,sizes,periodic,overlap,helper.getCommunicator());
>
> const auto &grid_view = grid.leafGridView();
> if (grid_view.comm().rank() == 0){
> std::cout << "[INFO]\tDomain size = " << grid.domainSize() << std::endl;
> std::cout << "[INFO]\t#cells = " << grid_view.size(0) << std::endl;
> }
>
> // make a local operator
> using LOP = Slopes<GridView>;
> LOP lop(grid_view);
>
> // make a grid operator
> using Con = Dune::PDELab::NoConstraints;
> using VBE = Dune::PDELab::istl::VectorBackend<>;
> using MBE = Dune::PDELab::istl::BCRSMatrixBackend<>;
> using FEM = Dune::PDELab::OPBLocalFiniteElementMap<double,double,0,dim,Dune::GeometryType::cube>;
> using GFS = Dune::PDELab::GridFunctionSpace<GridView,FEM,Con,VBE>;
> using GO = Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,double,double,double>;
> using U = typename GO::Traits::Domain;
>
> Con con;
> FEM fem;
> GFS gfs(grid_view,fem,con);
> MBE mbe(1);
> GO go(gfs,gfs,lop,mbe);
>
> // apply the grid operator
> U x(gfs,0.0);
> U r(gfs,0.0);
> go.residual(x,r);
>
> return 0;
> }
> catch (Dune::Exception &e){
> std::cerr << "Dune reported error: " << e << std::endl;
> }
> catch (...){
> std::cerr << "Unknown exception thrown!" << std::endl;
> }
> }
>
> _______________________________________________
> dune-pdelab mailing list
> dune-pdelab at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune-pdelab
--
Prof. Dr. Christian Engwer
Institut für Numerische und Angewandte Mathematik
Fachbereich Mathematik und Informatik der Universität Münster
Einsteinstrasse 62
48149 Münster
E-Mail christian.engwer at uni-muenster.de
Telefon +49 251 83-35067
FAX +49 251 83-32729
More information about the dune-pdelab
mailing list