[dune-pdelab] Finite Volume Method with slope reconstruction / alpha_volume_post_skeleton
Gregor Corbin
corbin at mathematik.uni-kl.de
Tue May 9 16:37:51 CEST 2017
Hi,
I misunderstood then, how the assembly works. Now that I added the
doSkeletonTwoSided flag, it seems to work.
However, my solution is still not ideal, since every slope is computed,
stored and loaded twice.
@Markus : I think we have a similar problem, namely that we need
(distant) neighbor information, where the local operator provides only
local information. My current workaround is somewhat ugly and works only
for explicit time-stepping. It involves writing a limiter class that
calls the residual method of the above local operator in its "void
preStage(X &x)" method. Maybe you could do something similar? If you are
interested, I can give you the code.
Cheers,
Gregor
On 09.05.2017 14:35, Christian Engwer wrote:
> 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
>
More information about the dune-pdelab
mailing list